Ported C-style 'dos_tetra' to C++

In order to wrap the ATM routines by Python using TRIQS wrapping
tools it is necessary to modify the interface to 'dos_tetra3d'.
The major changes involved replacing direct NumPy arrays with
TRIQS arrays which can be converted to Python arrays using library
tools.

Also, some small changes were necessary to port the functions
from C99 complex numbers to C++ style.

CMakeList is added to automatize building of the ATM library.

c++/plovasp/atm/CMakeLists.txt

@@ -0,0 +1,9 @@
+# Linking and include info
+add_library(atm_c dos_tetra3d.hpp dos_tetra3d.cpp argsort.hpp argsort.cpp)
+set_target_properties(atm_c PROPERTIES LINKER_LANGUAGE CXX)
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/c++/plovasp/atm ${TRIQS_INCLUDE_ALL})
+
+#add_executable(test_atm test2py.cpp)
+#target_link_libraries(test_atm atm_c)
+
+add_subdirectory(test)

c++/plovasp/atm/argsort.cpp

@@ -18,7 +18,7 @@
 * TRIQS. If not, see <http://www.gnu.org/licenses/>.
 *
 *******************************************************************************/
-#include <stdlib.h>
+#include <cstdlib>
 
 int cmp(const void *a, const void *b)
 {

c++/plovasp/atm/dos_tetra3d.cpp

@@ -18,15 +18,18 @@
 * TRIQS. If not, see <http://www.gnu.org/licenses/>.
 *
 *******************************************************************************/
-#include <Python.h>
+#include <triqs/arrays.hpp>
+#include <iostream>
+#include <complex>
 
-#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
+#include "argsort.hpp"
+#include "dos_tetra3d.hpp"
 
-#include <numpy/arrayobject.h>
+//#define __TETRA_DEBUG
-#include <complex.h>
+#define __TETRA_ARRAY_VIEW
-//#include "tetra.h"
+
-#include "argsort.h"
+using triqs::arrays::array;
-#include "dos_tetra3d.h"
+using triqs::arrays::array_view;
 
 /***************************************************
 
@@ -35,6 +38,23 @@ Lambin et al., PRB 29, 6, 3430 (1984).
 
 ***************************************************/
 
+/// Main function
+//#ifdef __TETRA_ARRAY_VIEW
+//void tet_dos3d(double en, array_view<double, 1>& eigk,
+//                 array_view<long, 2>& itt, int ntet,
+//                 array<double, 2>& cti);
+//#else
+//void tet_dos3d(double en, array<double, 1>& eigk,
+//                 array<long, 2>& itt, int ntet,
+//                 array<double, 2>& cti);
+//#endif
+
+
+/// Internal functions
+int dos_corner_weights(double en, double *eigs, int *inds, double *ci);
+int dos_tet_weights(double en, double *eigs, int *inds, double *ct);
+int dos_reorder(double en, double *e, int *inds);
+
 static double F(double en, double e1, double e2, double e3, double e4);
 static double K2(double en, double e1, double e2, double e3);
 static double K1(double en, double e1, double e2);
@@ -43,180 +63,141 @@ static void fun_dos_case1(double en, double *eigs, double *ci);
 static void fun_dos_case2(double en, double *eigs, double *ci);
 static void fun_dos_case3(double en, double *eigs, double *ci);
 
-int dos_corner_weights(double en, double *eigs, int *inds, double *ci);
+static const int NUM_TET_CORNERS = 4;
-int dos_tet_weights(double en, double *eigs, int *inds, double *ct);
+static const std::complex<double> I(0.0, 1.0);
+static const double small = 2.5e-2, tol = 1e-8;
 
-const static int NUM_TET_CORNERS = 4;
+/*
-
+  Returns corner contributions to the DOS of a band
-static PyMethodDef c_tetraMethods[] = {
+*/
-  {"dos_weights_3d", tetra_DOS3D, METH_VARARGS,
+#ifdef __TETRA_ARRAY_VIEW
-   "C-implementation of the tetrahedron method for calculating DOS"},
+array_view<double, 2> dos_tetra_weights_3d(array_view<double, 1> eigk, double en, array_view<long, 2> itt)
-  {NULL, NULL, 0, NULL}
+#else
-};
+array<double, 2> dos_tetra_weights_3d(array<double, 1> eigk, double en, array<long, 2> itt)
-
+#endif
-PyMODINIT_FUNC
-initc_atm_dos(void)
 {
-  (void) Py_InitModule("c_atm_dos", c_tetraMethods);
+  int ntet; /// Number of tetrahedra
-  import_array();
-}
-
-//
-// Integration_Weights
-//
-static PyObject *
-tetra_DOS3D(PyObject *self, PyObject *args)
-{
-
-  PyArrayObject *py_eigk, *py_itt;  // Input Numpy arrays
-  PyArrayObject *py_cti;  // Output Numpy array
-//  npy_int64 *itt;  // C-pointer to the 'itt' array
-  npy_intp *dims;  // Dimensions of the output array 'rti'
-  double *cti; // C-pointer to 'eigk' and 'rti' arrays
-  double e;
-//  double et[4]; // Real part of eigenvalues
-//  int strd_itt[2], strd_eigk, strd_cti[2], it_cti; // Strides
-  int ntet; // Indices of sorted 'eigs'; number of tetrahedra
 // Auxiliary variables and loop indices
-  int nd;
 
-// Data-preparation part
+  if (first_dim(itt) != NUM_TET_CORNERS + 1)
-  if (!PyArg_ParseTuple(args, "O!dO!", 
-      &PyArray_Type, &py_eigk, &e, &PyArray_Type, &py_itt))
-    return NULL;
-
-// Sanity tests (the types are assumed to be checked in the python wrapper)
-//  nd = py_eigk->nd;
-  nd = PyArray_NDIM(py_eigk);
-  if (nd != 1)
   {
-    PyErr_SetString(PyExc_ValueError, "  Array 'eigk' must be 1D");
+      TRIQS_RUNTIME_ERROR << "  The first dimension of 'itt' must be equal to 5";
-    return NULL;
   }
 
-//  nd = py_itt->nd;
+  ntet = second_dim(itt);
-  nd = PyArray_NDIM(py_itt);
-  if (nd != 2)
-  {
-    PyErr_SetString(PyExc_ValueError, "  Array 'itt' must be 2D");
-    return NULL;
-  }
 
-//  nd = py_itt->dimensions[0];
+  array<double, 2> cti(NUM_TET_CORNERS, ntet); // Corner weights to be returned
-  dims = PyArray_DIMS(py_itt);
-  if (dims[0] != 5)
-  {
-    PyErr_SetString(PyExc_ValueError, 
-       "  The first dimension of 'itt' must be equal to 5");
-    return NULL;
-  }
 
-//  ntet = (int) py_itt->dimensions[1];
+//  tet_dos3d(e, eigk, itt, ntet, cti);
-  ntet = dims[1];
-
-//  eigk = (double *)py_eigk->data;
-//  strd_eigk = py_eigk->strides[0] / sizeof(npy_float64);
-
-//  itt = (npy_int64 *)py_itt->data;
-//  strd_itt[0] = py_itt->strides[0] / sizeof(npy_int64);
-//  strd_itt[1] = py_itt->strides[1] / sizeof(npy_int64);
-
-// Resulting array (the question is whether dims is copied or not?)
-  cti = (double *)malloc(NUM_TET_CORNERS * ntet * sizeof(double));
-  dims = (npy_intp *)malloc(2 * sizeof(npy_intp));
-
-  dims[0] = NUM_TET_CORNERS;
-  dims[1] = ntet;
-
-  py_cti = (PyArrayObject *)PyArray_SimpleNewFromData(2, dims, NPY_DOUBLE, cti);
-
-//  strd_cti[0] = py_cti->strides[0] / sizeof(double);
-//  strd_cti[1] = py_cti->strides[1] / sizeof(double);
-
-// Main part: now we can fill the 'cti' array and it will be returned
-// by 'py_cti' as a numpy array
 
 //
-// Main function
+// Main algorithm (transferred from 'tet_dos3d()')
 //
-//  tet_dos3d(e, eigk, strd_eigk, itt, ntet, strd_itt, cti, strd_cti);
-  tet_dos3d(e, py_eigk, py_itt, ntet, py_cti);
- 
-  return PyArray_Return(py_cti);
-}
-
-//void tet_dos3d(double en, double *eigk, int strd_eigk,
-//                 npy_int64 *itt, int ntet, int *strd_itt,
-//                 double *cti, int *strd_cti)
-void tet_dos3d(double en, PyArrayObject *py_eigk, 
-                 PyArrayObject *py_itt, int ntet,
-                 PyArrayObject *py_cti)
-{
   double eigs[4], ci[4];
+  
   int i, it, ik, inds[4], flag;
-// **** DEBUG
+#ifdef __TETRA_DEBUG
   double ct, ci_sum;
+#endif
 
 // Loop over tetrahedra (triangles)
   for (it = 0; it < ntet; it++)
   {
-//    it_cti = it * strd_cti[1];
-// Sort eigenvalues and obtain indices of the sorted array
-//   eigs: sorted eigenvalues
-//   inds: index map
     for (i = 1; i < 5; i++) 
     {
-//      ik = itt[i * strd_itt[0] + it * strd_itt[1]];
+      ik = itt(i, it);
-      ik = ((int *)PyArray_GETPTR2(py_itt, i, it))[0];
+      eigs[i - 1] = eigk(ik);
-//      eigs[i - 1] = eigk[ik * strd_eigk];
-      eigs[i - 1] = ((double *)PyArray_GETPTR1(py_eigk, ik))[0];
     }
 
-// corner weights for a single tetrahedron
+// Corner weights for a single tetrahedron
     dos_corner_weights(en, eigs, inds, ci);
+
+#ifdef __TETRA_DEBUG
     for(i = 0, ci_sum = 0.0; i < 4; i++)
       ci_sum += ci[i];
 
     flag = dos_tet_weights(en, eigs, inds, &ct);
     if(fabs(ct - ci_sum) > tol)
       {
-        printf("  *** Error in weights: it = %d, flag = %d, en = %lf", it, flag, en);
+        std::cout << "  *** Error in weights: it = " << it <<" flag = " << flag << ", en = " << en;
         for(i = 0; i < 4; i++)
-          printf(", e[%d] = %lf", i, eigs[i]);
+          std::cout << ", e[" << i << "] = " << eigs[i];
-        printf(",    c_diff = %le\n", fabs(ct - ci_sum));
+        std::cout << ",    c_diff = " << fabs(ct - ci_sum) << std::endl;
-        return;
+
+        TRIQS_RUNTIME_ERROR << "  Failed consistency check";
       }  
-//    printf("  it = %d, flag = %d", it, j);
+#endif
-//    for(i = 0; i < 4; i++)
-//      printf(", e[%d] = %lf", i, eigs[i]);
-//    printf(",    c_diff = %le\n", fabs(ct - ci_sum));
 
-//    if(j < 4)
+  for(i = 0; i < 4; i++)
-//    {
-//      printf("  flag = %d, e = %lf", j, en);
-//      for(i = 0; i < 4; i++)
-//        printf(", eigs[%d] = %lf", i, eigs[i]);
-//      printf("\n");
-//      printf("  ci = ");
-//      for(i = 0; i < 4; i++)
-//        printf(", %lf", ci[i]);
-//      printf("\n");
-//    }
-
-    for(i = 0; i < 4; i++)
     {
-//      j = inds[i] * strd_cti[0] + it_cti;
+      cti(inds[i], it) = ci[i];
-//      cti[j] = ci[i];
-      ((double *)PyArray_GETPTR2(py_cti, inds[i], it))[0] = ci[i];
     }
 
   }  // it = 1, ntet
+
+  return array_view<double,2>(cti);
 }
 
+//#ifdef __TETRA_ARRAY_VIEW
+//void tet_dos3d(double en, array_view<double, 1>& eigk,
+//                 array_view<long, 2>& itt, int ntet,
+//                 array<double, 2>& cti)
+//#else
+//void tet_dos3d(double en, array<double, 1>& eigk,
+//                 array<long, 2>& itt, int ntet,
+//                 array<double, 2>& cti)
+//#endif
+//{
+//  double eigs[4], ci[4];
+//  
+//  int i, it, ik, inds[4], flag;
+//#ifdef __TETRA_DEBUG
+//  double ct, ci_sum;
+//#endif
+//
+//// Loop over tetrahedra (triangles)
+//  for (it = 0; it < ntet; it++)
+//  {
+//    for (i = 1; i < 5; i++) 
+//    {
+//      ik = itt(i, it);
+//      eigs[i - 1] = eigk(ik);
+//    }
+//
+//// Corner weights for a single tetrahedron
+//    dos_corner_weights(en, eigs, inds, ci);
+//
+//#ifdef __TETRA_DEBUG
+//    for(i = 0, ci_sum = 0.0; i < 4; i++)
+//      ci_sum += ci[i];
+//
+//    flag = dos_tet_weights(en, eigs, inds, &ct);
+//    if(fabs(ct - ci_sum) > tol)
+//      {
+//        std::cout << "  *** Error in weights: it = " << it <<" flag = " << flag << ", en = " << en;
+//        for(i = 0; i < 4; i++)
+//          std::cout << ", e[" << i << "] = " << eigs[i];
+//        std::cout << ",    c_diff = " << fabs(ct - ci_sum) << std::endl;
+//        return;
+//      }  
+//#endif
+//
+//  for(i = 0; i < 4; i++)
+//    {
+//      cti(inds[i], it) = ci[i];
+//    }
+//
+//  }  // it = 1, ntet
+//}
+ 
+/// Corner contributions to DOS
 int dos_corner_weights(double en, double *eigs, int *inds,
                    double *ci)
 {
   int flag, i;
+// Sort eigenvalues and obtain indices of the sorted array
+//   eigs: sorted eigenvalues
+//   inds: index map
   flag = dos_reorder(en, eigs, inds);
 
   switch(flag)
@@ -251,11 +232,14 @@ int dos_corner_weights(double en, double *eigs, int *inds,
   return flag;
 }
 
+/// Total (tetrahedron) contribution to DOS.
+/// Here, it is calculated directly using an analytical formula.
+/// This is mainly needed for debugging.
 int dos_tet_weights(double en, double *eigs, int *inds,
                     double *ct)
 {
   double e1, e2, e3, e4;
-  double complex s;
+  std::complex<double> s;
   int flag;
   flag = dos_reorder(en, eigs, inds);
 
@@ -275,9 +259,9 @@ int dos_tet_weights(double en, double *eigs, int *inds,
       s = fmin(fabs(e1 - e2), fabs(e3 - e1));
       s = fmin(fabs(s), fabs(e4 - e1));
       s /= 100.0;
-      s = fmax(s, 1.0e-20) * I;
+      s = fmax(std::abs(s), 1.0e-20) * I;
 
-      *ct = 3.0 * creal((en - e1 + s) * (en - e1 + s) / ((e2 - e1 + s) * (e3 - e1 + s) * (e4 - e1 + s)));
+      *ct = 3.0 * std::real((en - e1 + s) * (en - e1 + s) / ((e2 - e1 + s) * (e3 - e1 + s) * (e4 - e1 + s)));
     }
  
     break;
@@ -294,9 +278,9 @@ int dos_tet_weights(double en, double *eigs, int *inds,
       s = fmin(fabs(s), fabs(e4 - e1));
       s = fmin(fabs(s), fabs(e4 - e2));
       s /= 100.0;
-      s = fmax(s, 1.0e-20) * I;
+      s = fmax(std::abs(s), 1.0e-20) * I;
 
-      *ct = 3.0 * creal((
+      *ct = 3.0 * std::real((
           (e3 - en + s) * (en - e2 + s) / ((e4 - e2 + s) * (e3 - e2 + s) * (e3 - e1 + s)) +
           (e4 - en + s) * (en - e1 + s) / ((e4 - e1 + s) * (e4 - e2 + s) * (e3 - e1 + s))));
     }
@@ -311,9 +295,9 @@ int dos_tet_weights(double en, double *eigs, int *inds,
       s = fmin(fabs(e4 - e2), fabs(e4 - e1));
       s = fmin(fabs(s), fabs(e4 - e3));
       s /= 100.0;
-      s = fmax(s, 1.0e-20) * I;
+      s = fmax(std::abs(s), 1.0e-20) * I;
 
-      *ct = 3.0 * creal((e4 - en + s) * (e4 - en + s) / ((e4 - e1 + s) * (e4 - e2 + s) * (e4 - e3 + s)));
+      *ct = 3.0 * std::real((e4 - en + s) * (e4 - en + s) / ((e4 - e1 + s) * (e4 - e2 + s) * (e4 - e3 + s)));
     }
 
     break;
@@ -333,6 +317,8 @@ int dos_tet_weights(double en, double *eigs, int *inds,
   return flag;
 }
 
+/// Sorts eigenvalues and also determines eigenvalue degeneracies.
+/// Returns a case number corresponding to a combination of degeneracies.
 int dos_reorder(double en, double *e, int *inds)
 {
   double *ptrs[4], e_tmp[4];
@@ -358,13 +344,6 @@ int dos_reorder(double en, double *e, int *inds)
 static void fun_dos_case1(double en, double *eigs, double *ci)
 {
   double e1, e2, e3, e4;
-//  int i;
-
-//  if(fabs(eigs[1] - eigs[0]) < tol)
-//  {
-//    for(i = 0; i < 4; i++) ci[i] = 0.0;
-//    return;
-//  }
 
   e1 = eigs[0];
   e2 = eigs[1];
@@ -386,12 +365,6 @@ static void fun_dos_case2(double en, double *eigs, double *ci)
 {
   double e1, e2, e3, e4;
 
-//  if(fabs(eigs[2] - eigs[1]) < tol)
-//  {
-//    for(i = 0; i < 4; i++) ci[i] = 0.0;
-//    return;
-//  }
-
   e1 = eigs[0];
   e2 = eigs[1];
   e3 = eigs[2];
@@ -438,12 +411,6 @@ static void fun_dos_case3(double en, double *eigs, double *ci)
 {
   double e1, e2, e3, e4;
 
-//  if(fabs(eigs[3] - eigs[2]) < tol)
-//  {
-//    for(i = 0; i < 4; i++) ci[i] = 0.0;
-//    return;
-//  }
-
   e1 = eigs[0];
   e2 = eigs[1];
   e3 = eigs[2];
@@ -463,49 +430,53 @@ static void fun_dos_case3(double en, double *eigs, double *ci)
 
 static double F(double en, double e1, double e2, double e3, double e4)
 {
-  double complex s;
+  std::complex<double> s;
 
   if(fabs(e1 - e3) > tol && fabs(e4 - e2) > tol)
     return (e1 - en) * (en - e2) / ((e1 - e3) * (e4 - e2));
   else
   {
+// Regularization to avoid division by zero
     s = fmin(fabs(e3 - e1), fabs(e4 - e2));
     s /= 100.0;
-    s = fmax(s, 1.0e-20) * I;
+    s = fmax(std::abs(s), 1.0e-20) * I;
    
-    return creal((e1 - en + s) * (en - e2 + s) / ((e1 - e3 + s) * (e4 - e2 + s)));
+    return std::real((e1 - en + s) * (en - e2 + s) / ((e1 - e3 + s) * (e4 - e2 + s)));
   }
 }
 
 static double K2(double en, double e1, double e2, double e3)
 {
-  double complex s;
+  std::complex<double> s;
 
   if(fabs(e1 - e3) > tol && fabs(e1 - e2) > tol)
     return (en - e1) / ((e2 - e1) * (e3 - e1));
   else
   {
+// Regularization to avoid division by zero
     s = fmin(fabs(e3 - e1), fabs(e1 - e2));
     s /= 100.0;
-    s = fmax(s, 1.0e-20) * I;
+    s = fmax(std::abs(s), 1.0e-20) * I;
    
-    return creal((en - e1 + s) / ((e2 - e1 + s) * (e3 - e1 + s)));
+    return std::real((en - e1 + s) / ((e2 - e1 + s) * (e3 - e1 + s)));
   }
 }
 
 static double K1(double en, double e1, double e2)
 {
-  double complex s;
+  std::complex<double> s;
 
   if(fabs(e1 - e2) > tol)
     return (e1 - en) / ((e2 - e1) * (e2 - e1));
   else
   {
+// Regularization to avoid division by zero
     s = fabs(e1 - e2);
     s /= 100.0;
-    s = fmax(s, 1.0e-20) * I;
+    s = fmax(std::abs(s), 1.0e-20) * I;
    
-    return creal((e1 - en + s) / ((e2 - e1 + s) * (e2 - e1 + s)));
+    return std::real((e1 - en + s) / ((e2 - e1 + s) * (e2 - e1 + s)));
   }
 }
 
+

c++/plovasp/atm/dos_tetra3d.h

@@ -1,42 +0,0 @@
-/*******************************************************************************
-*
-* This file is part of the ATM library.
-*
-* Copyright (C) 2010 by O. E. Peil
-*
-* TRIQS is free software: you can redistribute it and/or modify it under the
-* terms of the GNU General Public License as published by the Free Software
-* Foundation, either version 3 of the License, or (at your option) any later
-* version.
-*
-* TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
-* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
-* details.
-*
-* You should have received a copy of the GNU General Public License along with
-* TRIQS. If not, see <http://www.gnu.org/licenses/>.
-*
-*******************************************************************************/
-#ifndef __C_DOS_TETRA3D_H__
-#define __C_DOS_TETRA3D_H__
-
-#include <Python.h>
-#include <numpy/arrayobject.h>
-
-static PyObject *tetra_DOS3D(PyObject *self, PyObject *args);
-
-//void tet_dos3d(double en, double *eigk, int strd_eigk,
-//                 npy_int64 *itt, int ntet, int *strd_itt,
-//                 double *cti, int *strd_cti);
-void tet_dos3d(double en, PyArrayObject *py_eigk,
-                 PyArrayObject *py_itt, int ntet,
-                 PyArrayObject *py_cti);
-int dos_corner_weights(double en, double *eigs, int *inds,
-                 double *ci);
-
-int dos_reorder(double en, double *e, int *inds);
-
-static const double small = 2.5e-2, tol = 1e-8;
-
-#endif

c++/plovasp/atm/dos_tetra3d.hpp

@@ -0,0 +1,46 @@
+/*******************************************************************************
+*
+* This file is part of the ATM library.
+*
+* Copyright (C) 2010 by O. E. Peil
+*
+* TRIQS is free software: you can redistribute it and/or modify it under the
+* terms of the GNU General Public License as published by the Free Software
+* Foundation, either version 3 of the License, or (at your option) any later
+* version.
+*
+* TRIQS is distributed in the hope that it will be useful, but WITHOUT ANY
+* WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
+* FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
+* details.
+*
+* You should have received a copy of the GNU General Public License along with
+* TRIQS. If not, see <http://www.gnu.org/licenses/>.
+*
+*******************************************************************************/
+#pragma once
+
+#include <triqs/arrays.hpp>
+
+using triqs::arrays::array;
+using triqs::arrays::array_view;
+
+/// DOS of a band by analytical tetrahedron method
+///
+///   Returns corner weights for all tetrahedra for a given band and real energy.
+array_view<double, 2> 
+dos_tetra_weights_3d(array_view<double, 1> eigk, /// Band energies for each k-point
+                     double en, /// Energy at which DOS weights are to be calculated
+                     array_view<long, 2> itt /// Tetrahedra defined by k-point indices
+);
+//array<double, 2> 
+//dos_tetra_weights_3d(array<double, 1> eigk, /// Band energies for each k-point
+//                     double e, /// Energy at which DOS weights are to be calculated
+//                     array<long, 2> itt /// Tetrahedra defined by k-point indices
+//);
+int dos_corner_weights(double en, double *eigs, int *inds, double *ci);
+int dos_tet_weights(double en, double *eigs, int *inds, double *ct);
+int dos_reorder(double en, double *e, int *inds);
+
+
+