The core test is converting numpy array into matrix 5x5 in Armadillo C++ inversion matrix into matrix C and verification myMat * C = E , afterwards we can safely return matrix C to Python Runtime Module
I was able to succeed with armadillo-10.6.1 setup on fedora 34 only via build from source , previously installed
$ sudo dnf install cmake openblas-devel lapack-devel \
arpack-devel SuperLU-devel
$ tar -xvf armadillo-10.6.1.tar.xz
$ cd armadillo-10.6.1
$ ./configure
$ make
$ sudo make install
Here we follow [1] all updates were colored in blue inside code
(.env) [boris@fedora33server NUMPYCPP]$ cat examplePlus.cpp
#define NPY_NO_DEPRECATED_API NPY_1_9_API_VERSION
extern "C" {
#include <Python.h>
#include <numpy/arrayobject.h>
}
#include <exception>
#include <cassert>
#include <string>
#include <type_traits>
#include <map>
#include <vector>
#include <armadillo>
class WrongDimensions : public std::exception
{
public:
WrongDimensions() {}
const char* what() const noexcept { return msg.c_str(); }
private:
std::string msg = "The dimensions were incorrect";
};
class NotImplemented : public std::exception
{
public:
NotImplemented() {}
const char* what() const noexcept { return msg.c_str(); }
private:
std::string msg = "Not implemented";
};
class BadArrayLayout : public std::exception
{
public:
BadArrayLayout() {}
const char* what() const noexcept { return msg.c_str(); }
private:
std::string msg = "The matrix was not contiguous";
};
static const std::vector<npy_intp> getPyArrayDimensions(PyArrayObject* pyarr)
{
npy_intp ndims = PyArray_NDIM(pyarr);
npy_intp* dims = PyArray_SHAPE(pyarr);
std::vector<npy_intp> result;
for (int i = 0; i < ndims; i++) {
result.push_back(dims[i]);
}
return result;
}
/* Checks the dimensions of the given array. Pass -1 for either dimension to say you don't
* care what the size is in that dimension. Pass dimensions (X, 1) for a vector.
*/
static bool checkPyArrayDimensions(PyArrayObject* pyarr, const npy_intp dim0, const npy_intp dim1)
{
const auto dims = getPyArrayDimensions(pyarr);
assert(dims.size() <= 2 && dims.size() > 0);
if (dims.size() == 1) {
return (dims[0] == dim0 || dim0 == -1) && (dim1 == 1 || dim1 == -1);
}
else {
return (dims[0] == dim0 || dim0 == -1) && (dims[1] == dim1 || dim1 == -1);
}
}
template<typename outT>
static arma::Mat<outT> convertPyArrayToArma(PyArrayObject* pyarr, int nrows, int ncols)
{
if (!checkPyArrayDimensions(pyarr, nrows, ncols)) throw WrongDimensions();
int arrTypeCode;
if (std::is_same<outT, uint16_t>::value) {
arrTypeCode = NPY_UINT16;
}
else if (std::is_same<outT, double>::value) {
arrTypeCode = NPY_DOUBLE;
}
else {
throw NotImplemented();
}
const auto dims = getPyArrayDimensions(pyarr);
if (dims.size() == 1) {
outT* dataPtr = static_cast<outT*>(PyArray_DATA(pyarr));
return arma::Col<outT>(dataPtr, dims[0], true);
}
else {
PyArray_Descr* reqDescr = PyArray_DescrFromType(arrTypeCode);
if (reqDescr == NULL) throw std::bad_alloc();
PyArrayObject* cleanArr = (PyArrayObject*)PyArray_FromArray(pyarr, reqDescr, NPY_ARRAY_FARRAY);
if (cleanArr == NULL) throw std::bad_alloc();
reqDescr = NULL;
outT* dataPtr = static_cast<outT*>(PyArray_DATA(cleanArr));
arma::Mat<outT> result (dataPtr, dims[0], dims[1], true);
Py_DECREF(cleanArr);
return result;
}
}
static PyObject* convertArmaToPyArray(const arma::mat& matrix)
{
npy_intp ndim = matrix.is_colvec() ? 1 : 2;
npy_intp nRows = static_cast<npy_intp>(matrix.n_rows); // NOTE: This narrows the integer
npy_intp nCols = static_cast<npy_intp>(matrix.n_cols);
npy_intp dims[2] = {nRows, nCols};
PyObject* result = PyArray_SimpleNew(ndim, dims, NPY_DOUBLE);
if (result == NULL) throw std::bad_alloc();
double* resultDataPtr = static_cast<double*>(PyArray_DATA((PyArrayObject*)result));
for (int i = 0; i < nRows; i++) {
for (int j = 0; j < nCols; j++) {
resultDataPtr[i * nCols + j] = matrix(i, j);
}
}
return result;
}
extern "C" {
static PyObject* example_testFunction(PyObject* self, PyObject* args)
{
int flag = 1;
PyArrayObject* myArray = NULL;
if (!PyArg_ParseTuple(args, "iO!", &flag, &PyArray_Type, &myArray)) {
return NULL;
}
PyObject* output = NULL;
try {
arma::mat myMat = convertPyArrayToArma<double>(myArray, -1, -1);
arma::mat A = {{ 1, 2, 3, 2, 4},
{ 2, 3, 4, 1, 5},
{ 3, 4, 5, 7, 1},
{ 4, 5, 6, 1, 2},
{ 5, 6, 7, 3, 1}};
arma::mat myOut;
if ( flag == 1) {
myOut = myMat * A ;
}
if ( flag == 2) {
myOut = arma::inv(myMat);
}
if (flag == 0) {
arma::mat C= arma::inv(myMat);
myOut = myMat * C ;
}
output = convertArmaToPyArray(myOut);
}
PyErr_NoMemory();
Py_XDECREF(output);
return NULL;
}
catch (const std::exception& err) {
PyErr_SetString(PyExc_RuntimeError, err.what());
Py_XDECREF(output);
return NULL;
}
return output;
}
static PyMethodDef example_methods[] =
{
{"test_function", example_testFunction, METH_VARARGS, "A test function"},
{NULL, NULL, 0, NULL}
};
static struct PyModuleDef example_module = {
PyModuleDef_HEAD_INIT,
"example", /* name of module */
NULL, /* module documentation, may be NULL */
-1, /* size of per-interpreter state of the module,
or -1 if the module keeps state in global variables. */
example_methods
};
PyMODINIT_FUNC
PyInit_example(void)
{
import_array();
PyObject* m = PyModule_Create(&example_module);
if (m == NULL) return NULL;
return m;
}
}
