[petsc-users] petsc4py: PetSc is no longer initialised when calling C Function from Cython

[Jeff Wiens]

The point of the example was to perform the PetSc operation in C.
Although you could do the dot product in Cython or Python, it defeats
my purpose. I am attempting to write a wrapper for an existing PetSc C
code base. I would like to setup PetSc and construct the initial
conditions from python. I would then like to send this information to
the PetSc C program to execute and then return to Python the final
state as a petsc4py vector. If there is a better way of doing this, I
would love to hear it. However, because of the size of the existing
PetSc C code, rewriting it in Cython or Python is not an option.

Jeff

On Thu, Oct 20, 2011 at 1:26 PM, Ethan Coon <ecoon at lanl.gov> wrote:
> On Thu, 2011-10-20 at 13:05 -0700, Jeff Wiens wrote:
>> Eventually, I would like
>> to Initialise PETSc and construct my PETSc vectors in python and then
>> have my C code perform operations on them. Again, this yet another
>> test problem.
>>
>
> Maybe I don't completely understand your planned working model, but I
> don't think it's the right choice.
>
> petsc4py is very good at making setup and "high level" features of PETSc
> work well quickly and easily.  If you plan to call things like
> VecSetSizes() from C, you might as well write your entire program in C.
> It would be much easier to simply call VecSetSizes() from petsc4py  via
> PETSc.Vec().setSizes().
>
> What I suspect you're trying to avoid is doing the work in python,
> including things like numerical calculations on the data contained
> within the Vec.  Doing this in C via Cython is a good idea.  The better
> model for implementing this sort of thing is to do all of your setup in
> python, getting the data from the Vec, and passing that data off to your
> cython/c code, i.e. something like:
>
> -- Python Code --
>
> import sys, petsc4py
> petsc4py.init(sys.argv)
> from petsc4py import PETSc
>
> v1 = PETSc.Vec().createSeq(3)
> v2 = v1.duplicate()
> v1.setFromOptions()
> v2.setFromOptions()
>
> v1.set(3.)
> v2.set(6.)
>
> import dot
> val = dot.mydot(v1[...], v2[...]) # note that these ellipses treat the underlying C-data stored in Vec as data to a numpy array
>
> v1.destroy()
> v2.destroy()
>
>
>
> -- Cython code --
> import numpy
> cimport numpy
>
> DTYPE = numpy.float
> ctypedef numpy.float_t DTYPE_t
>
> cdef mydot( numpy.ndarray[DTYPE_t, ndim=1] v1,
>            numpy.ndarray[DTYPE_t, ndim=1] v2):
>    cdef int np = v1.shape[0]
>    cdef lcv
>    cdef DTYPE_t val
>
>    val = 0.0
>    for lcv in range(np):
>        val += v1[lcv]*v2[lcv]
>    return val
>
> ---------
>
> Ethan
>
>
>
>
>> The code is as follows:
>>
>> ------- Python Code -----------
>>
>> import sys, petsc4py
>> petsc4py.init(sys.argv)
>> from petsc4py import PETSc
>>
>> import dot
>> val = dot.mydot()
>> PETSc.Sys.Print( "PETSc Dot Product: %i"%val )
>>
>> ------- Cython Code -----------
>>
>> cdef extern from "c_files/dot.h":
>>    void init()
>>    double dot()
>>
>> def mydot():
>>    #init();
>>    return dot();
>>
>> ------- C Code -----------
>>
>> void init()
>> {
>>     PetscInitialize(NULL,NULL,(char *)0,NULL);
>> }
>>
>>
>> double dot()
>> {
>>     PetscScalar r;
>>     Vec v1,v2;
>>
>>     VecCreate(PETSC_COMM_WORLD, &v1);
>>     VecCreate(PETSC_COMM_WORLD, &v2);
>>
>>     VecSetSizes(v1, PETSC_DECIDE, 5000);
>>     VecSetSizes(v2, PETSC_DECIDE, 5000);
>>     VecSetFromOptions(v1);
>>     VecSetFromOptions(v2);
>>
>>     VecSet(v1,1.0);
>>     VecSet(v2,2.0);
>>
>>     VecAssemblyBegin(v1);
>>     VecAssemblyEnd(v1);
>>     VecAssemblyBegin(v2);
>>     VecAssemblyEnd(v2);
>>
>>     VecDot(v1, v2, &r);
>>     return r;
>> }
>
> --
> ------------------------------------
> Ethan Coon
> Post-Doctoral Researcher
> Applied Mathematics - T-5
> Los Alamos National Laboratory
> 505-665-8289
>
> http://www.ldeo.columbia.edu/~ecoon/
> ------------------------------------
>
>