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

[Ethan Coon]

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/
------------------------------------