3

I am writing a Fortran code and I would like to use some special functions and methods from Python libraries. This is a Python code:

from mpmath import *
from scipy.optimize import *

def g(A,B,t):   
   return newton(lambda x: (x - A*polylog(1.5, B*exp(-t*x))), 0.0)

In fortran code I would like to pass real values A,B,t and get in return the value of g(A,B,t) (possibly complex). This can also be done in C.

I would like to mention that mixing Python with other languages is something new to me.

Bociek
  • 1,195
  • 2
  • 13
  • 28
  • 2
    You should probably start here [Embedding Python in Another Application (python docs)](https://docs.python.org/3.5/extending/embedding.html). – Thomas Orozco Dec 18 '15 at 11:59
  • Always use tag [tag:fortran] and only add the version when necessary to distinguish that your question is specific. For example that you cannot use Fortran 2008 but only Fortran 90. – Vladimir F Героям слава Dec 18 '15 at 12:34
  • See http://stackoverflow.com/questions/17198642/callback-python-from-fortran http://stackoverflow.com/questions/21836670/how-to-expose-python-callbacks-to-fortran-using-modules http://stackoverflow.com/questions/20265697/how-to-return-a-value-from-a-python-callback-in-fortran-using-f2py – Vladimir F Героям слава Dec 18 '15 at 12:36
  • 3
    Do you realize that a lot of SciPy and related libraries is implemented in Fortran and C already? Would you consider just calling Fortran or C from Fortran? Do you have experience calling C from Fortran? If so you can use the Python C API. – John Zwinck Dec 18 '15 at 15:28
  • @JohnZwinck - that is my goal to use Python C API and then write a simple Fortran wrapper. – Bociek Dec 18 '15 at 16:12
  • 1
    You might want to consider to yourself and mark that post as the correct answer ;) – Mr_Pouet Dec 19 '15 at 00:03

1 Answers1

1

SOLUTION

In case someone is interested I have solved this problem. The following was of great help: stackoverflow and yolinux.

rootC.c

#include "rootC.h"
#include <Python.h>
#include <stdlib.h>

void Initialize ()
{
    Py_Initialize();
}

void Finalize ()
{
    Py_Finalize();
}

void root_(double* A, double* B, double* t, double* x)
{
    PyObject *pName, *pModule, *pFunc;
    PyObject *pArgs, *pValue, *sys, *path;
    static int i;

    if (i == 0)
    {    
        ++i;
        Initialize();
        atexit(Finalize);
    }

    sys  = PyImport_ImportModule("sys");
    path = PyObject_GetAttrString(sys, "path");
    PyList_Append(path, PyString_FromString("."));

    pName = PyString_FromString("rootPY");
    pModule = PyImport_Import(pName);

    if (!pModule)
    {
        PyErr_Print();
        printf("ERROR in pModule\n");
        exit(1);
    }

    pFunc = PyObject_GetAttrString(pModule, "root");
    pArgs = PyTuple_New(3);
    PyTuple_SetItem(pArgs, 0, PyFloat_FromDouble((*A)));
    PyTuple_SetItem(pArgs, 1, PyFloat_FromDouble((*B)));
    PyTuple_SetItem(pArgs, 2, PyFloat_FromDouble((*t)));
    
    pValue = PyObject_CallObject(pFunc, pArgs);
    *x     = PyFloat_AsDouble(pValue);
}

rootC.h

#ifndef ROOT_H_
#define ROOT_H_

void Initialize();
void Finalize();
void root_(double*, double*, double*, double*);

#endif

rootPY.py

from mpmath import polylog, exp
from scipy.optimize import newton

def root(A,B,t):   
    return abs(newton(lambda x: (x - A*polylog(1.5, B*exp(-t*x))), 0.0))

rootF.F90

program main
   implicit none
   real*8 A, B, t, x
   A = 0.4d0
   B = 0.3d0
   t = 0.1d0

   call root(A,B,t,x)

   write(*,*) x

end program main

main.c

#include "rootC.h"
#include <stdio.h>

int main() 
{
    double A = 0.4, B = 0.3, t = 0.1, x = 0.0;

    root_(&A,&B,&t,&x);

    printf("A = %f, B = %f, t = %f\n", A, B, t);

    printf("x = %.15f\n", x);

    return 0;
}

Makefile

CC = gcc
FC = gfortran
CFLAGS = -I/usr/include/python2.7
LFLAGS = -L/usr/local/lib -lpython2.7 -lm

.PHONY: all clean

all: rootF main

rootF: rootF.o rootC.o
    $(FC) $^ -o $@ $(LFLAGS)

rootF.o: rootF.F90
    $(FC) -c $< -o $@

main: main.o rootC.o
    $(CC) $^ -o $@ $(LFLAGS)

main.o: main.c
    $(CC) $(CFLAGS) -c $< -o $@

rootC.o: rootC.c
    $(CC) $(CFLAGS) -c $< -o $@

clean:
    rm -f *.o
Community
  • 1
  • 1
Bociek
  • 1,195
  • 2
  • 13
  • 28
  • Unfortunatrly it gives me Segmentation Fault when I call the function twice - at pModule = Py_Import... – Bociek Dec 19 '15 at 12:41
  • Solved thanks to [this](http://stackoverflow.com/questions/34371221/segmentation-fault-when-calling-c-function-with-python-c-api-twice/34375287#34375287) – Bociek Dec 19 '15 at 22:04
  • The only thing which need to be taken care is the reference count of Python objects. – Bociek Dec 20 '15 at 11:03