[petsc-dev] PetscSF in Fortran

[Barry Smith]

> On Oct 20, 2017, at 12:31 PM, Jed Brown <jed at jedbrown.org> wrote:
> 
> Barry Smith <bsmith at mcs.anl.gov> writes:
> 
>>   Adrian,
>> 
>>    You should not use F90Array1d *rptr as arguments in the Fortran interface. You should just use a regular Fortran one dimensional array of real/scalar.
>> Fortran doesn't handle polymorphism in this way at all. You have to have multiple f90 interfaces, one for each type and provide a C stub for real, int, or whatever else you want to send.
> 
> Barry, look at the "use mpi_f08" way of calling MPI from Fortran.

  Jed,

   Rather terse response. 

   Are you suggesting in PETSc we use type(*) to manage multiple types through the same function?  Looks doable, I wasn't aware of this. This could possibly reduce a lot of code duplication we currently have. 

   Still I would like to get rid of the use PetscDataType rather than write new code that uses it.

   I need to think more in this case.

   Waiting to hear from Adrian what types he needs to pass around (use of PetscDataType restricts to built in MPI datatypes regardless of what Fortran interface approach we use


  Barry

> 
>>   What types do you need to send? It is probably easier if we just write the Fortran interfaces for you.
>> 
>>   If you can send a simple Fortran PETS  code that uses PetscSFBcastBegin() and PetscSFBcastEnd() that I can use for testing then I will write them. Do a bcast with int and another with real/scalar.
>> 
>> 
>>   Barry
>> 
>> 
>> 
>>> On Oct 19, 2017, at 8:59 PM, Adrian Croucher <a.croucher at auckland.ac.nz> wrote:
>>> 
>>> hi
>>> 
>>> On 19/10/17 06:45, Matthew Knepley wrote:
>>>> On Tue, Oct 17, 2017 at 11:35 PM, Adrian Croucher <a.croucher at auckland.ac.nz> wrote:
>>>> 
>>>> 
>>>> So, now I'm trying to add Fortran bindings for PetscSFBcastBegin() and PetscSFBcastEnd().
>>>> 
>>>> From the C side I have added the following into src/vec/is/sf/interface/f90-custom/zsff90.c:
>>>> 
>>>> PETSC_EXTERN void PETSC_STDCALL petscsfbcastbegin_(PetscSF *sf, MPI_Datatype *unit, F90Array1d *rptr, F90Array1d *lptr , int *ierr PETSC_F90_2PTR_PROTO(rptrd) PETSC_F90_2PTR_PROTO(lptrd))
>>>> {
>>>>  PetscDataType ptype;
>>>>  const void* rootdata;
>>>>  void* leafdata;
>>>> 
>>>>  *ierr = PetscMPIDataTypeToPetscDataType(*unit, &ptype);if (*ierr) return;
>>>>  *ierr = F90Array1dAccess(rptr, ptype, (void**) &rootdata PETSC_F90_2PTR_PARAM(rptrd));if (*ierr) return;
>>>>  *ierr = F90Array1dAccess(lptr, ptype, (void**) &leafdata PETSC_F90_2PTR_PARAM(lptrd));if (*ierr) return;
>>>> 
>>>>  *ierr = PetscSFBcastBegin(*sf, *unit, rootdata, leafdata);
>>>> 
>>>> }
>>>> 
>>>> and similarly for petscsfbcastend_(). Does this look plausible?
>>>> 
>>>> Then some wrappers need to be added to src/vec/f90-mod/petscis.h90. I am not sure how to do those.
>>>> 
>>>> The difficulty is in declaring the arrays that are passed in, which can be of various types. In C they are declared as void*, but I'm not sure what to do with that in Fortran. I can't seem to find any other example wrappers in PETSc to model it on either. Any suggestions?
>>> 
>>> 
>>> I think this is working now by just declaring those void* C variables as type(*) in the Fortran interface, e.g.:
>>> 
>>>      Interface
>>>         Subroutine PetscSFBcastBegin(sf,unit,rarray,
>>>     &       larray,ierr)
>>>           use petscisdef
>>>           PetscSF :: sf
>>>           PetscInt :: unit
>>>           type(*) :: rarray(:)
>>>           type(*) :: larray(:)
>>>           PetscErrorCode :: ierr
>>>         End Subroutine PetscSFBcastBegin
>>>      End Interface
>>> 
>>> The only difficulty I have left with this is in the MPI_Datatype variable. I'd forgotten that these datatypes are all different in C and Fortran as well.
>>> 
>>> I amended the C interface code to the following, to convert the Fortran MPI datatype (actually an integer) to a C MPI_Datatype:
>>> 
>>> PETSC_EXTERN void PETSC_STDCALL petscsfbcastbegin_(PetscSF *sf, MPI_Fint *unit, F90Array1d *rptr, F90Array1d *lptr , int *ierr PETSC_F90_2PTR_PROTO(rptrd) PETSC_F90_2PTR_PROTO(lptrd))
>>> {
>>>  PetscDataType pdtype;
>>>  MPI_Datatype dtype;
>>>  const void* rootdata;
>>>  void* leafdata;
>>> 
>>>  dtype = MPI_Type_f2c(*unit);
>>>  *ierr = PetscMPIDataTypeToPetscDataType(dtype, &pdtype);if (*ierr) return;
>>>  *ierr = F90Array1dAccess(rptr, pdtype, (void**) &rootdata PETSC_F90_2PTR_PARAM(rptrd));if (*ierr) return;
>>>  *ierr = F90Array1dAccess(lptr, pdtype, (void**) &leafdata PETSC_F90_2PTR_PARAM(lptrd));if (*ierr) return;
>>> 
>>>  *ierr = PetscSFBcastBegin(*sf, dtype, rootdata, leafdata);
>>> 
>>> }
>>> 
>>> The problem is this only seems to work if I declare the datatype in the calling Fortran code to be of the appropriate C MPI datatype, e.g. MPI_INT, rather than the corresponding Fortran datatype, e.g. MPI_INTEGER (which causes PetscMPIDataTypeToPetscDataType() to fail, as something weird gets passed in for dtype).
>>> 
>>> I was expecting the opposite to be true. It doesn't seem right to have to use the C datatypes in Fortran code (confusing if the Fortran datatypes are used elsewhere). So I suspect I've messed something up. Anyone have any ideas?
>>> 
>>> - Adrian
>>> -- 
>>> Dr Adrian Croucher
>>> Senior Research Fellow
>>> Department of Engineering Science
>>> University of Auckland, New Zealand
>>> email: 
>>> a.croucher at auckland.ac.nz
>>> 
>>> tel: +64 (0)9 923 4611
>>>