Skip to content
Snippets Groups Projects
Commit 9cb49a8c authored by Jonathan Lambrechts's avatar Jonathan Lambrechts
Browse files

python : fullMatrix<int> <=> numpy int array

parent 868a9817
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Numeric/fullMatrix.cpp
+ 55
− 0
View file @ 9cb49a8c
...@@ -66,6 +66,16 @@ void fullVector<std::complex<double> >::setAll(const fullVector<std::complex<dou ...@@ -66,6 +66,16 @@ void fullVector<std::complex<double> >::setAll(const fullVector<std::complex<dou
F77NAME(zcopy)(&_r, m._data, &stride, _data, &stride); F77NAME(zcopy)(&_r, m._data, &stride, _data, &stride);
} }
template<>
void fullMatrix<int>::setAll(const fullMatrix<int> &m)
{
if (_r != m._r || _c != m._c )
Msg::Fatal("fullMatrix size does not match");
int N = _r * _c;
for (int i = 0; i < N; ++i)
_data[i] = m._data[i];
}
template<> template<>
void fullMatrix<double>::setAll(const fullMatrix<double> &m) void fullMatrix<double>::setAll(const fullMatrix<double> &m)
{ {
...@@ -213,6 +223,13 @@ void fullMatrix<double>::multWithATranspose(const fullVector<double> &x, double ...@@ -213,6 +223,13 @@ void fullMatrix<double>::multWithATranspose(const fullVector<double> &x, double
&beta, y._data, &INCY); &beta, y._data, &INCY);
} }
template<>
void fullMatrix<int>::scale(const double s)
{
for(int i = 0; i < _r * _c; ++i) _data[i] *= s;
}
#endif #endif
...@@ -585,3 +602,41 @@ bool fullMatrix<double>::invert(fullMatrix<double> &result) const ...@@ -585,3 +602,41 @@ bool fullMatrix<double>::invert(fullMatrix<double> &result) const
} }
#endif #endif
template<>
void fullMatrix<double>::print(const std::string name, const std::string format) const
{
std::string rformat = (format == "") ? "%12.5E " : format;
int ni = size1();
int nj = size2();
printf("double %s [ %d ][ %d ]= { \n", name.c_str(),ni,nj);
for(int I = 0; I < ni; I++){
printf("{ ");
for(int J = 0; J < nj; J++){
printf(rformat.c_str(), (*this)(I, J));
if (J!=nj-1)printf(",");
}
if (I!=ni-1)printf("},\n");
else printf("}\n");
}
printf("};\n");
}
template<>
void fullMatrix<int>::print(const std::string name, const std::string format) const
{
std::string rformat = (format == "") ? "%12d " : format;
int ni = size1();
int nj = size2();
printf("int %s [ %d ][ %d ]= { \n", name.c_str(),ni,nj);
for(int I = 0; I < ni; I++){
printf("{ ");
for(int J = 0; J < nj; J++){
printf(rformat.c_str(), (*this)(I, J));
if (J!=nj-1)printf(",");
}
if (I!=ni-1)printf("},\n");
else printf("}\n");
}
printf("};\n");
}
Numeric/fullMatrix.h
+ 3
− 18
View file @ 9cb49a8c
...@@ -430,7 +430,7 @@ class fullMatrix ...@@ -430,7 +430,7 @@ class fullMatrix
if (init0) if (init0)
setAll(scalar(0.)); setAll(scalar(0.));
} }
fullMatrix(int r, int c, double *data) fullMatrix(int r, int c, scalar *data)
: _r(r), _c(c), _data(data), _own_data(false) : _r(r), _c(c), _data(data), _own_data(false)
{ {
setAll(scalar(0.)); setAll(scalar(0.));
...@@ -526,7 +526,7 @@ class fullMatrix ...@@ -526,7 +526,7 @@ class fullMatrix
_own_data = false; _own_data = false;
_data = original._data + c_start * _r; _data = original._data + c_start * _r;
} }
void setAsProxy(double *data, int r, int c) void setAsProxy(scalar *data, int r, int c)
{ {
if(_data && _own_data) if(_data && _own_data)
delete [] _data; delete [] _data;
...@@ -774,22 +774,7 @@ class fullMatrix ...@@ -774,22 +774,7 @@ class fullMatrix
#endif #endif
; ;
void print(const std::string name = "", const std::string format = "%12.5E ") const void print(const std::string name = "", const std::string format = "") const;
{
int ni = size1();
int nj = size2();
printf("double %s [ %d ][ %d ]= { \n", name.c_str(),ni,nj);
for(int I = 0; I < ni; I++){
printf("{ ");
for(int J = 0; J < nj; J++){
printf(format.c_str(), (*this)(I, J));
if (J!=nj-1)printf(",");
}
if (I!=ni-1)printf("},\n");
else printf("}\n");
}
printf("};\n");
}
void binarySave (FILE *f) const void binarySave (FILE *f) const
{ {
......
wrappers/gmshpy/gmshNumeric.i
+ 16
− 0
View file @ 9cb49a8c
...@@ -28,10 +28,26 @@ ...@@ -28,10 +28,26 @@
%include "JacobianBasis.h" %include "JacobianBasis.h"
%ignore fullMatrix<double>::operator()(int, int); %ignore fullMatrix<double>::operator()(int, int);
%ignore fullVector<double>::operator()(int); %ignore fullVector<double>::operator()(int);
%ignore fullMatrix<int>::operator()(int, int);
%ignore fullMatrix<int>::eig;
%ignore fullMatrix<int>::luSubstitute;
%ignore fullMatrix<int>::invert;
%ignore fullMatrix<int>::invertInPlace;
%ignore fullMatrix<int>::axpy;
%ignore fullMatrix<int>::gemm;
%ignore fullMatrix<int>::mult;
%ignore fullMatrix<int>::luSolve;
%ignore fullMatrix<int>::luFactor;
%ignore fullMatrix<int>::svd;
%ignore fullMatrix<int>::multAddy;
%ignore fullMatrix<int>::multWithATranspose;
%ignore fullMatrix<int>::multOnBlock;
%ignore fullMatrix<int>::determinant;
%include "fullMatrix.h" %include "fullMatrix.h"
%include "simpleFunction.h" %include "simpleFunction.h"
%template(fullMatrixDouble) fullMatrix<double>; %template(fullMatrixDouble) fullMatrix<double>;
%template(fullVectorDouble) fullVector<double>; %template(fullVectorDouble) fullVector<double>;
%template(fullMatrixInt) fullMatrix<int>;
%template(simpleFunctionDouble) simpleFunction<double>; %template(simpleFunctionDouble) simpleFunction<double>;
%include "simpleFunctionPython.h" %include "simpleFunctionPython.h"
%include "nodalBasis.h" %include "nodalBasis.h"
......
wrappers/gmshpy/gmshtypemaps.i
+ 212
− 3
View file @ 9cb49a8c
...@@ -61,7 +61,7 @@ ...@@ -61,7 +61,7 @@
return fm; return fm;
%#ifdef HAVE_NUMPY %#ifdef HAVE_NUMPY
PyArrayObject *array = (PyArrayObject*)obj; PyArrayObject *array = (PyArrayObject*)obj;
if (PyArray_Check(array) && PyArray_ISFARRAY_RO(array) && PyArray_NDIM(array) == 2) { if (PyArray_Check(array) && PyArray_ISFARRAY_RO(array) && PyArray_NDIM(array) == 2 && PyArray_DESCR(array)->kind == 'f') {
newMatrix = true; newMatrix = true;
return new fullMatrix<double>((double*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1)); return new fullMatrix<double>((double*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1));
} }
...@@ -75,7 +75,7 @@ ...@@ -75,7 +75,7 @@
newMatrix = true; newMatrix = true;
return fm; return fm;
} }
fullMatrix<double> *objToFullMatrixRW(PyObject *obj, bool &newMatrix) fullMatrix<double> *objToFullMatrixRW(PyObject *obj, bool &newMatrix)
{ {
fullMatrix<double> *fm = NULL; fullMatrix<double> *fm = NULL;
...@@ -84,7 +84,7 @@ ...@@ -84,7 +84,7 @@
return fm; return fm;
%#ifdef HAVE_NUMPY %#ifdef HAVE_NUMPY
PyArrayObject *array = (PyArrayObject*)obj; PyArrayObject *array = (PyArrayObject*)obj;
if (PyArray_Check(array) && PyArray_ISFARRAY(array) && PyArray_NDIM(array) == 2) { if (PyArray_Check(array) && PyArray_ISFARRAY(array) && PyArray_NDIM(array) == 2 && PyArray_DESCR(array)->kind == 'f') {
newMatrix = true; newMatrix = true;
return new fullMatrix<double>((double*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1)); return new fullMatrix<double>((double*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1));
} }
...@@ -152,6 +152,154 @@ ...@@ -152,6 +152,154 @@
} }
%#endif %#endif
} }
%fragment("fullMatrixConversionInt", "header", fragment="fullMatrixConversionInit") {
%#undef HAVE_DLOPEN
%#include "fullMatrix.h"
%#include "GmshConfig.h"
%#ifdef HAVE_NUMPY
%#define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
%#include <numpy/arrayobject.h>
%#if NPY_API_VERSION < 0x00000007
%#define NPY_ARRAY_FARRAY NPY_FARRAY
%#define NPY_ARRAY_ALIGNED NPY_ALIGNED
%#define NPY_ARRAY_F_CONTIGUOUS NPY_F_CONTIGUOUS
%#define NPY_ARRAY_WRITEABLE NPY_WRITEABLE
%#define PyArray_SetBaseObject(a, b) PyArray_BASE(a) = b
%#endif
%#endif
fullMatrix<int> *pySequenceToFullMatrixInt(PyObject *o) {
fullMatrix<int> *fm = NULL;
if (!PySequence_Check(o))
return NULL;
long int nRow = PySequence_Length(o);
for (int i = 0; i < PySequence_Length(o); ++i) {
PyObject *l = PySequence_GetItem(o, i);
if (!PySequence_Check(l)){
if (fm != NULL)
delete fm;
return NULL;
}
long int nCol = PySequence_Length(l);
if (i == 0)
fm = new fullMatrix<int>(nRow, nCol);
else if (nCol != fm->size2()) {
delete fm;
return NULL;
}
for (int j = 0; j < nCol; ++j) {
PyObject *v = PySequence_GetItem(l, j);
if (PyLong_Check(v)) {
(*fm)(i, j) = (int) PyLong_AsLong(v);
}
else {
delete fm;
return NULL;
}
}
}
return fm;
}
fullMatrix<int> *objToFullMatrixIntRO(PyObject *obj, bool &newMatrix, PyObject *&tmpObject)
{
fullMatrix<int> *fm = NULL;
SWIG_ConvertPtr(obj,(void **) &fm, SWIGTYPE_p_fullMatrixT_int_t, 1);
if (fm)
return fm;
%#ifdef HAVE_NUMPY
PyArrayObject *array = (PyArrayObject*)obj;
if (PyArray_Check(array) && PyArray_ISFARRAY_RO(array) && PyArray_NDIM(array) == 2 && PyArray_DESCR(array)->kind == 'i') {
newMatrix = true;
return new fullMatrix<int>((int*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1));
}
if ((tmpObject = PyArray_FROMANY(obj, NPY_INT, 2, 2, NPY_ARRAY_FARRAY))) {
array = (PyArrayObject*)tmpObject;
newMatrix = true;
return new fullMatrix<int>((int*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1));
}
%#endif
if ((fm = pySequenceToFullMatrixInt(obj)))
newMatrix = true;
return fm;
}
fullMatrix<int> *objToFullMatrixIntRW(PyObject *obj, bool &newMatrix)
{
fullMatrix<int> *fm = NULL;
SWIG_ConvertPtr(obj,(void **) &fm, SWIGTYPE_p_fullMatrixT_int_t, 1);
if (fm)
return fm;
%#ifdef HAVE_NUMPY
PyArrayObject *array = (PyArrayObject*)obj;
if (PyArray_Check(array) && PyArray_ISFARRAY(array) && PyArray_NDIM(array) == 2 && PyArray_DESCR(array)->kind == 'i') {
newMatrix = true;
return new fullMatrix<int>((int*)PyArray_DATA(array), PyArray_DIM(array, 0), PyArray_DIM(array, 1));
}
%#endif
return NULL;
}
%#ifdef HAVE_NUMPY
%#if PY_MAJOR_VERSION==2 && PY_MINOR_VERSION==6
static void deleteCapsuleArrayInt(void *capsule)
{
delete [](double*)(capsule);
}
%#else
static void deleteCapsuleArrayInt(PyObject *capsule)
{
delete [](double*)PyCapsule_GetPointer(capsule, NULL);
}
%#endif
PyObject *fullMatrix2PyArray(fullMatrix<int> &fm)
{
npy_intp dims[2] = {fm.size1(), fm.size2()};
int *data = &fm.operator()(0, 0);
/*PyObject *array = PyArray_New(&PyArray_Type, 2, dims, NPY_DOUBLE, NULL, NULL, 0, NPY_ARRAY_F_CONTIGUOUS, NULL);
// copy data
memcpy((void*)PyArray_DATA(array), data, dims[0] * dims[1] * sizeof(int));*/
// do not copy data
PyObject *array = PyArray_New(&PyArray_Type, 2, dims, NPY_INT, NULL, (void*)data, 0, NPY_ARRAY_F_CONTIGUOUS, NULL);
PyArray_UpdateFlags((PyArrayObject*)array, NPY_ARRAY_ALIGNED);
if (fm.getOwnData()) {
fm.setOwnData(false);
%#if PY_MAJOR_VERSION==2 && PY_MINOR_VERSION==6
PyObject *capsule = PyCObject_FromVoidPtr((void*)data, deleteCapsuleArrayInt);
%#else
PyObject *capsule = PyCapsule_New((void*) data, NULL, deleteCapsuleArrayInt);
%#endif
PyArray_SetBaseObject((PyArrayObject*)array, capsule);
}
return array;
}
PyObject *fullMatrix2PyArrayConst(const fullMatrix<int> &fm)
{
npy_intp dims[2] = {fm.size1(), fm.size2()};
int *data = (int*)fm.getDataPtr();
// do not copy data
PyObject *array = PyArray_New(&PyArray_Type, 2, dims, NPY_INT, NULL, (void*)data, 0, NPY_ARRAY_F_CONTIGUOUS, NULL);
PyArray_UpdateFlags((PyArrayObject*)array, NPY_ARRAY_ALIGNED);
return array;
}
PyObject *fullMatrix2PyArrayProxy(const fullMatrix<int> &fm)
{
npy_intp dims[2] = {fm.size1(), fm.size2()};
int *data = (int*)fm.getDataPtr();
PyObject *array = PyArray_New(&PyArray_Type, 2, dims, NPY_INT, NULL, (void*)data, 0, NPY_ARRAY_F_CONTIGUOUS, NULL);
PyArray_UpdateFlags((PyArrayObject*)array, NPY_ARRAY_ALIGNED |NPY_ARRAY_WRITEABLE);
return array;
}
PyObject *fullVector2PyArrayProxy(fullVector<int> &fv)
{
npy_intp dims[1] = {fv.size()};
int *data = &fv.operator()(0);
PyObject *array = PyArray_New(&PyArray_Type, 1, dims, NPY_INT, NULL, (void*)data, 0, NPY_ARRAY_F_CONTIGUOUS, NULL);
PyArray_UpdateFlags((PyArrayObject*)array, NPY_ARRAY_ALIGNED | NPY_ARRAY_WRITEABLE);
return array;
}
%#endif
}
%typemap(in, fragment="fullMatrixConversion") const fullMatrix<double> &(PyObject *tmpObject = NULL, bool newMatrix = false){ %typemap(in, fragment="fullMatrixConversion") const fullMatrix<double> &(PyObject *tmpObject = NULL, bool newMatrix = false){
$1 = objToFullMatrixRO($input, newMatrix, tmpObject); $1 = objToFullMatrixRO($input, newMatrix, tmpObject);
...@@ -165,6 +313,19 @@ ...@@ -165,6 +313,19 @@
if (newMatrix$argnum && $1) delete $1; if (newMatrix$argnum && $1) delete $1;
} }
%typemap(in, fragment="fullMatrixConversionInt") const fullMatrix<int> &(PyObject *tmpObject = NULL, bool newMatrix = false){
$1 = objToFullMatrixIntRO($input, newMatrix, tmpObject);
if (!$1) {
PyErr_Format(PyExc_TypeError, "cannot convert argument %i to a fullMatrix<int>", $argnum);
SWIG_fail;
}
}
%typemap(freearg) const fullMatrix<int> &{
if (tmpObject$argnum) Py_DECREF(tmpObject$argnum);
if (newMatrix$argnum && $1) delete $1;
}
%typemap(in, fragment="fullMatrixConversion") fullMatrix<double> &(bool newMatrix = false){ %typemap(in, fragment="fullMatrixConversion") fullMatrix<double> &(bool newMatrix = false){
$1 = objToFullMatrixRW($input, newMatrix); $1 = objToFullMatrixRW($input, newMatrix);
if (!$1) { if (!$1) {
...@@ -176,6 +337,17 @@ ...@@ -176,6 +337,17 @@
if (newMatrix$argnum && $1) delete $1; if (newMatrix$argnum && $1) delete $1;
} }
%typemap(in, fragment="fullMatrixConversionInt") fullMatrix<int> &(bool newMatrix = false){
$1 = objToFullMatrixIntRW($input, newMatrix);
if (!$1) {
PyErr_Format(PyExc_TypeError, "cannot convert argument %i to a writable fullMatrix<int>", $argnum);
SWIG_fail;
}
}
%typemap(freearg) fullMatrix<int> &{
if (newMatrix$argnum && $1) delete $1;
}
%typemap(in, fragment="fullMatrixConversion") const fullMatrix<double> *(PyObject *tmpObject = NULL, bool newMatrix = false){ %typemap(in, fragment="fullMatrixConversion") const fullMatrix<double> *(PyObject *tmpObject = NULL, bool newMatrix = false){
$1 = objToFullMatrixRO($input, newMatrix, tmpObject); $1 = objToFullMatrixRO($input, newMatrix, tmpObject);
if (!$1) { if (!$1) {
...@@ -188,6 +360,18 @@ ...@@ -188,6 +360,18 @@
if (newMatrix$argnum && $1) delete $1; if (newMatrix$argnum && $1) delete $1;
} }
%typemap(in, fragment="fullMatrixConversionInt") const fullMatrix<int> *(PyObject *tmpObject = NULL, bool newMatrix = false){
$1 = objToFullMatrixIntRO($input, newMatrix, tmpObject);
if (!$1) {
PyErr_Format(PyExc_TypeError, "cannot convert argument %i to a fullMatrix<int>", $argnum);
SWIG_fail;
}
}
%typemap(freearg) const fullMatrix<int> *{
if (tmpObject$argnum) Py_DECREF(tmpObject$argnum);
if (newMatrix$argnum && $1) delete $1;
}
%typemap(in, fragment="fullMatrixConversion") fullMatrix<double> *(bool newMatrix = false){ %typemap(in, fragment="fullMatrixConversion") fullMatrix<double> *(bool newMatrix = false){
$1 = objToFullMatrixRW($input, newMatrix); $1 = objToFullMatrixRW($input, newMatrix);
if (!$1) { if (!$1) {
...@@ -199,6 +383,17 @@ ...@@ -199,6 +383,17 @@
if (newMatrix$argnum && $1) delete $1; if (newMatrix$argnum && $1) delete $1;
} }
%typemap(in, fragment="fullMatrixConversionInt") fullMatrix<int> *(bool newMatrix = false){
$1 = objToFullMatrixIntRW($input, newMatrix);
if (!$1) {
PyErr_Format(PyExc_TypeError, "cannot convert argument %i to a writable fullMatrix<int>", $argnum);
SWIG_fail;
}
}
%typemap(freearg) fullMatrix<int> *{
if (newMatrix$argnum && $1) delete $1;
}
#ifdef HAVE_NUMPY #ifdef HAVE_NUMPY
%typemap(out, fragment="fullMatrixConversion") fullMatrix<double> { %typemap(out, fragment="fullMatrixConversion") fullMatrix<double> {
$result = fullMatrix2PyArray($1); $result = fullMatrix2PyArray($1);
...@@ -215,5 +410,19 @@ ...@@ -215,5 +410,19 @@
%typemap(out, fragment="fullMatrixConversion") fullVector<double>& { %typemap(out, fragment="fullMatrixConversion") fullVector<double>& {
$result = fullVector2PyArrayProxy(*$1); $result = fullVector2PyArrayProxy(*$1);
} }
%typemap(out, fragment="fullMatrixConversionInt") fullMatrix<int> {
$result = fullMatrix2PyArray($1);
}
%typemap(out, fragment="fullMatrixConversionInt") const fullMatrix<int>& {
$result = fullMatrix2PyArrayConst(*$1);
}
%typemap(out, fragment="fullMatrixConversionInt") fullMatrix<int>& {
$result = fullMatrix2PyArrayProxy(*$1);
}
#endif #endif
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment