Using scipy in Fortran

Using scipy in Fortran

In this example, we will solve an optimization problem with the help of scipy in Fortran. We want to solve the following linear problem (from https://de.wikipedia.org/wiki/Lineare_Optimierung, in German):

minimize f(x1, x2) = -300 * x1 -500 * x2

x1 + 2*x2 <= 170
x1 +   x2 <= 150
     3*x2 <= 180

x1 >= 0
x2 >= 0

We will use the function scipy.optimize.linprog from the scipy package (assuming it is installed). In Python the code would be something like this:

import numpy as np
import scipy.optimize as opt

c = np.array([-300., -500.])
b_ub = np.array([170., 150., 180.])

A_ub = np.array([[1., 2.],
                 [1., 1.],
                 [0., 3.]])

retval = opt.linprog(c, A_ub, b_ub)
x = retval.x
objective_fun_value = retval.fun

print("Solution: x = {0}".format(x))
print("Value of objective function: fun = {0}".format(
       objective_fun_value))

Below we find the roughly equivalent code in Fortran with forpy (without error handling for now). We have to import scipy.optimize, create numpy wrappers for 3 Fortran arrays, build the argument tuple for linprog and then call the function. linprog returns its results as an object with various attributes. We use the getattribute method on the object to retrieve a result that interest us. Then we cast the retrieved attribute to a Fortran value.

program scipy_example01
use forpy_mod
implicit none

integer, parameter:: dp = kind(0.d0)

integer :: ierror
type(module_py) :: opt
type(ndarray) :: nd_c, nd_b_ub, nd_A_ub, nd_x
type(object) :: retval, attr
type(tuple) :: args

real(dp) :: c(2) = [-300._dp, -500._dp]
real(dp) :: b_ub(3) = [170._dp, 150._dp, 180._dp]

real(dp) :: A_ub(3,2)
real(dp), dimension(:), pointer :: x
real(dp) :: objective_fun_value

A_ub(1,1) = 1.0_dp 
A_ub(2,1) = 1.0_dp 
A_ub(3,1) = 0.0_dp 

A_ub(1,2) = 2.0_dp
A_ub(2,2) = 1.0_dp
A_ub(3,2) = 3.0_dp

ierror = forpy_initialize()
ierror = import_py(opt, "scipy.optimize")

ierror = ndarray_create(nd_c, c)
ierror = ndarray_create(nd_b_ub, b_ub)
ierror = ndarray_create(nd_A_ub, A_ub)

ierror = tuple_create(args, 3)
ierror = args%setitem(0, nd_c)
ierror = args%setitem(1, nd_A_ub)
ierror = args%setitem(2, nd_b_ub)

ierror = call_py(retval, opt, "linprog", args)

ierror = retval%getattribute(attr, "x")
ierror = cast(nd_x, attr)
ierror = nd_x%get_data(x)
call attr%destroy

ierror = retval%getattribute(attr, "fun")
ierror = cast(objective_fun_value, attr)
call attr%destroy

write(*,*) "Solution: x = ", x
write(*,*) "Value of objective function: fun = ", objective_fun_value

call retval%destroy
call args%destroy
call nd_c%destroy
call nd_b_ub%destroy
call nd_A_ub%destroy
call nd_x%destroy
call opt%destroy

call forpy_finalize()

end program