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.

6

1 Answer 1

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
Sign up to request clarification or add additional context in comments.

3 Comments

Unfortunatrly it gives me Segmentation Fault when I call the function twice - at pModule = Py_Import...
Solved thanks to this
The only thing which need to be taken care is the reference count of Python objects.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.