Re: [pyamg-user] Re: Calling pyamg from an external Fortran90 application
Ondrej Certik
>
>
> On Mon, Oct 5, 2009 at 11:54 AM, Ondrej Certik <ond...@certik.cz> wrote:
>>
>> On Mon, Oct 5, 2009 at 5:03 AM, Luke <luke....@gmail.com> wrote:
>> >
>> > AS,
>> >
>> > Thanks for you comments. Glad you've liked the package so far.
>> >
>> > As for your Fortran code, we're working on a pure C-branch
>> > specifically for this reason. It's a ways off, however. I've tried
>> > both of your solutions:
>> >
>> > Solution1:
>> > For me this worked very well. I was impressed with F2Py. I had to do
>> > very little coding to interface Python with Fortran and did not notice
>> > a performance drop. This is an easy first pass that's worth trying.
>> >
>> > Solution 2:
>> > I've never been able to do this :) Embedding the interpreter works
>> > (for me) for very small python examples, but nothing with SciPy/NumPy/
>> > PyAMG. You may be a little more savvy, however; I didn't spend too
>> > much time investigating.
>>
>> The example how to do embed python (and numpy+scipy works!) into a C
>> (resp. C++) app is here:
>>
>> http://hpfem.org/hermes1d/
>>
>> e.g. do:
>>
>> git clone http://hpfem.math.unr.edu/git/hermes1d.git
>> cd hermes1d
>> python -c 'print "set(WITH_PYTHON yes)\n"' > CMake.vars
>> cmake .
>> make
>> cd examples/schroedinger
>> ./schroedinger
>>
>>
>> and a matplotlib window will popup, I use numpy for solving the
>> eigenproblem an this is all done from *within* a C++ app, see the
>> main.cpp (see the line 86). Here is an example of how it looks like if
>> you run it:
>>
>>
>> http://groups.google.com/group/hermes1d/browse_thread/thread/3f73e5bde7ebe527
>>
>> the hermes2d is a very young code, so you will almost surely discover
>> some bugs, but the python <-> C interface works great and without any
>> problems, I use it in other projects too.
>>
>> Ondrej
>>
>
> Ondrej,
>
> Nice example. One difficulty, however, is that the data is passed by
> writing to a file and then rereading in the python script. In the other
That is just a quick hack. You can ignore that.
> example, a full copy is performed. The real challenge for PyAMG would be to
> pass a reference to the data as data to a numpy object. If you have any
> example of this, that would be helpful.
It's super easy to pass a reference to the data, *without* copying
anything. See below how to do it.
In the other example, this Cython code is responsible for converting
C++ array into a numpy array:
cdef ndarray array_double_c2numpy(double *A, int len):
from numpy import empty
cdef ndarray vec = empty([len], dtype="double")
cdef double *pvec = <double *>vec.data
memcpy(pvec, A, len*sizeof(double))
return vec
(see python/hermes1d/_hermes1d.pyx, line 94), so currently the data is
copied, because that is the safest approach (you never know if the
C++/Python code doesn't free the memory if the other one needs it). In
my case the memcopy is really fast and not the bottleneck.
If however, you know that you can safely use the exact same memory
from both C++ and Python, or the memcpy is really killing you, just
use this patch:
diff --git a/python/hermes1d/_hermes1d.pyx b/python/hermes1d/_hermes1d.pyx
index c1d8cdd..3586069 100644
--- a/python/hermes1d/_hermes1d.pyx
+++ b/python/hermes1d/_hermes1d.pyx
@@ -94,8 +94,9 @@ cdef ndarray array_int_c2numpy(int *A, int len):
cdef ndarray array_double_c2numpy(double *A, int len):
from numpy import empty
cdef ndarray vec = empty([len], dtype="double")
- cdef double *pvec = <double *>vec.data
- memcpy(pvec, A, len*sizeof(double))
+ vec.data = <char *>A
+ #cdef double *pvec = <double *>vec.data
+ #memcpy(pvec, A, len*sizeof(double))
return vec
This doesn't copy any memory and I just tested that it works for me.
To get data from Python into C++ without copying, just use this code
(see the line 101) to convert from python to C++:
cdef api void array_double_numpy2c_inplace(object A_n, double **A_c, int *n):
cdef ndarray A = A_n
if not (A.nd == 1 and A.strides[0] == sizeof(double)):
from numpy import array
A = array(A.flat, dtype="double")
n[0] = len(A)
A_c[0] = <double *>(A.data)
you can see that no memory is copied, except if the numpy array is not
contiguous. I tested that and it works for me.
Note that this inplace creating of numpy arrays doesn't usually work
for me, since I use a code like this:
insert_int("len", len);
double _mat[len*len];
for(int i=0; i<len; i++)
for(int j=0; j<len; j++)
_mat[i*len+j] = mat->get(i, j);
insert_double_array("mat", _mat, len*len);
the insert_double_array() calls the array_double_c2numpy() above and
as you can see, the _mat memory gets destroyed after the
insert_double_array(), so one has to be really careful and think about
this all the time. To fix the above, I would have to do:
diff --git a/examples/schroedinger/main.cpp b/examples/schroedinger/main.cpp
index 50c35bf..6bfd107 100644
--- a/examples/schroedinger/main.cpp
+++ b/examples/schroedinger/main.cpp
@@ -37,7 +37,7 @@ double rhs(int num, double *x, double *weights,
void insert_matrix(DenseMatrix *mat, int len)
{
insert_int("len", len);
- double _mat[len*len];
+ double *_mat = new double[len*len];
for(int i=0; i<len; i++)
for(int j=0; j<len; j++)
_mat[i*len+j] = mat->get(i, j);
then it works. But as I said, it's a major pain to handle this memory by hand.
Ondrej
P.S. CCing cython-dev, maybe someone there could offer a better advise.
Ondrej Certik
One problem that I can see with this is that the "vec.data" is first
allocated, and then I should free it myself in the above function, and
then numpy will free the new "vec.data". This also takes care of:
> then it works. But as I said, it's a major pain to handle this memory by hand.
^^^ numpy actually frees it for me I think.
The remaining problem is how to make numpy not to allocate the memory
at all. But in any case, that imho is a minor problem, it should not
slow things down and if you free it yourself in the above function, it
should get you started.
Please ping me if you figure out how to force numpy not to allocate
the "vec.data" in the first place.
Ondrej
Ondrej Certik
Ok, I have fixed that now, the solution is to use the numpy C/API. So
if you do now:
git clone http://hpfem.math.unr.edu/git/hermes1d.git
cd hermes1d
python -c 'print "set(WITH_PYTHON yes)\n"' > CMake.vars
cmake .
make
cd examples/schroedinger
./schroedinger
it will be copying the data using memcpy. If you apply this patch:
diff --git a/examples/schroedinger/main.cpp b/examples/schroedinger/main.cpp
index 4800440..5d3be6e 100644
--- a/examples/schroedinger/main.cpp
+++ b/examples/schroedinger/main.cpp
@@ -36,12 +36,12 @@ double rhs(int num, double *x, double *weights,
void insert_matrix(DenseMatrix *mat, int len)
{
- double _mat[len*len];
+ double *_mat = new double[len*len];
for(int i=0; i<len; i++)
for(int j=0; j<len; j++)
_mat[i*len+j] = mat->get(i, j);
insert_object("len", c2py_int(len));
- insert_object("mat", c2numpy_double(_mat, len*len));
+ insert_object("mat", c2numpy_double_inplace(_mat, len*len));
cmd("_ = mat.reshape((len, len))");
cmd("del len");
cmd("del mat");
it will not be copying anything, nor allocating any memory in vain
during the conversion.
Ondrej