Investigate IoC

[aprokop]

The Inversion of Control in Futility is organized as follows.

User create some subroutine, for example

    SUBROUTINE testcalcResid() BIND(C)
      !Do some matvec
    ENDSUBROUTINE testcalcResid

This function is passed through a proxy by using C_FUNLOC:

call JFNK_Init(C_FUNLOC(testcalcResid), ...)

JFNK_Init has the following prototype:

    subroutine JFNK_Init(fptr,...) bind(C,NAME="JFNK_Init")
      IMPORT :: C_FUNPTR
      TYPE(C_FUNPTR),INTENT(IN),VALUE :: fptr
      ...
    endsubroutine

with a corresponding C interfaces:

extern "C" void JFNK_Init(int &id, void (*funptr)(), const int idx,
                          const int idF) {
    id = jfnk->new_data(funptr, evec->get_vec(idx), evec->get_vec(idF));
}

eventually ending in

class ModelEvaluator_JFNK : public NOX::Epetra::Interface::Required {
public:
ModelEvaluator_JFNK(void (*functionptr)(),Teuchos::RCP<Epetra_Vector> x, Teuchos::RCP<Epetra_Vector> F){
    fptr=functionptr;
    xloc=x;
    Floc=F;
}

...

bool computeF(const Epetra_Vector& x,
              Epetra_Vector& f,
              NOX::Epetra::Interface::Required::FillType)
    {
      // Residual calculation
      *xloc = x;
      fptr();
      f = *Floc;
      return true;
    }
private:
    void (*fptr)() = NULL;
    ...
};

[aprokop]

Here are the next steps to be done (building up on work in #14):

• [ ] Implement this pattern for operator in TrilinosHandle
• [ ] Implement this pattern for preconditioner in TrilinosHandle

It may be of use to CAM that does not build operator and preconditioner explicitly.

[aprokop]

I've written a prototype (available here). I've got it to compile, and call the Fortran function from C++, but it segfaults. It seems that on the Fortran side it gets messed up parameters. I wonder if I first need to bind the IoC function to some C function.

[aprokop]

@youngmit @bscollin It would be super nice if you could help me out with this.

[aprokop]

Fixed it, works now. Note to self: when passing double * from C/C++ side to Fortran, the corresponding variable should be

  real(C_double), dimension(:), value(out) :: x(*)

Notice the (*) at the end, this makes all the difference, and means that the variable may have arbitrary length.

[aprokop]

Couple approaches discussed with @tjfulle

1. Fixed function name approach (Abaqus umat)

   A user write a function with a well defined signature. This signature could include some void* data and some name. So it really could be multiple functions in one with where the exact function is chosen by name. In Trilinos linear solvers, for example, this could mean that a user could write a single "apply" call with a simple signature like fortran_apply(double*, double*, int n, char* name) and then this function would be called from C++. This is kind of similar to what ACME does for nonlinear solvers, where it provides calc_f, update_jac_state, etc, and does bind(C) to those functions on the Fortran side.

   The main benefit of this approach is the lack of need to pass around the function pointers, as the name of the function on the Fortran side is fixed. Instead, a user passes a (optional?) name and data. The likely drawback is that C++ will need to provide a dummy implementation of said function, which would then have to be replaced by a Fortran user provided function during the linking stage.

2. Function pointer approach

   A Fortran user creates a function and then passes a pointer to that function to the C++ side. The function signature is probably going to be exactly the same as in

   The benefits of this approach seem the lack of linking issues of the previous approach, and the fact that a user is not required to put all functions in the same place.

@kevans32 @sethrj I'd like to hear your thoughts on these 2 approaches.

What data to put in function API?

An orthogonal question to the whole discussion is what the function signature should be. In simple terms, should we strive to provide a signature that takes in ForTrilinos data (TpetraMultiVector, for example), or would Fortran users prefer the native data? If former, how would one implement this? @kevans32 what would ACME's preferred approach would be if we were able to provide, say, Jacobian which would operate on TpetraMultiVector?

I think the concrete steps we should start taking are:

• [ ] Create a small example which does the single function approach with POD data, and demonstrate the possibility of compilation when a) a user provides the implementation, and b) a user does not provide implementation of the target function. What are the challenges, and workarounds to make this work?
• [ ] Create a small example that does the function pointer approach with POD data.
• [ ] Create a small example using either approach which would implement a function on Fortran side with a signature that takes in ForTrilinos data like TpetraMultiVector.

@tjfulle What do you think?

[sethrj]

Lol I just realized most of the above stuff was written in 2017 not just recently. Here's what I had written. But I think in regard to the IOC, with SWIG we can generate a class that uses Fortran inheritance to allow the user to override and translate data, using the "director" capability for cross-language polymorphism. So for the C++ class NOX::Epetra::Interface::Required, we would have SWIG code like

%module(directors="1") example
%{
#include "NOX_Whatever_decl.hpp"
%}

%feature("director") NOX::Epetra::Interface::Required;

%include "NOX_Whatever_decl.hpp"

SWIG would generate a NOX::Epetra::Interface::Required class that includes conversion code for input and output of the variables, so that a Fortran daughter class can override those methods as if they were native fortran code:

type(jfnk_model), extends(Nox_Interface)  
  procedure :: computeF => my_compute_f
end type
contains
subroutine my_compute_f(x, f)
  type(Epetra_Vector), intent(in) :: x
  type(Epetra_Vector), intent(in,out) :: f
  integer(kind(FillType)) :: fill
  ! do things
end subroutine

SWIG would generate C++ code that implements the virtual override of the C++ computeF and would dispatch it to the Fortran overridable function. There would be some overhead from the extra function pointer calls and conversion involved, but I don't think it would be any worse than the best that could be done by processing it manually.

Cross referencing with https://github.com/sethrj/swig/issues/8

---

@aprokop Wish I'd been there for this discussion -- we could have gone over the SWIG "director" capability (which is implemented in multiple other languages, but currently not in Fortran). It basically allows you to implement a C++ virtual function call in the target language, with all the wrapping niceties in between.

Incidentally, did I show you the %bindc feature for C classes? From an input

extern "C" void JFNK_Init(int &id, void (*funptr)(), const int idx,
                          const int idF) {
    id = jfnk->new_data(funptr, evec->get_vec(idx), evec->get_vec(idF));
}

it will generate

subroutine JFNK_Init(fptr,...) bind(C,NAME="JFNK_Init")
      IMPORT :: C_FUNPTR
      TYPE(C_FUNPTR),INTENT(IN),VALUE :: fptr
      ...
end subroutine