homework 3: embedded fortran and C in Python

cfunc.c

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

void c_add(double* a, double* b, double* c)
{
	*c = *a + *b;
}

cfunc.h

void c_add(double* a, double* b, double* c);

gfunc.f90

module gfunc_module
implicit none
contains
subroutine multi(a, b, c)
    double precision, intent(in) :: a, b
    double precision, intent(out) :: c

    c = a * b
    
end subroutine multi
    
end module gfunc_module

merge.pyx

cdef extern from "cfunc.h":
	void c_add(double* a, double* b, double* c)

cdef extern from "pygfunc.h":
	void c_multi(double* a, double* b, double* c)

def pymulti(double a, double b):
	cdef:
		double c
	c_multi(&a, &b, &c)
	return c

def pyadd(double a, double b):
	cdef:
		double c
	c_add(&a, &b, &c)
	return c

pygfunc.f90

module gfunc_interface
use iso_c_binding, only: c_double
use gfunc_module, only: multi
implicit none
contains
subroutine c_multi(a, b, c) bind(c)
    real(c_double), intent(in) :: a, b
    real(c_double), intent(out) :: c
    call multi(a, b, c)
    
end subroutine c_multi
    
end module gfunc_interface

pygfunc.h

extern void c_multi(double* a, double* b, double* c);

setupm.py

from distutils.core import setup
from distutils.extension import Extension
from Cython.Distutils import build_ext
# This line only needed if building with NumPy in Cython file.
from numpy import get_include
from os import system

# compile the fortran modules without linking
fortran_mod_comp = 'gfortran gfunc.f90 -c -o gfunc.o -O3 -fPIC'
print(fortran_mod_comp)
system(fortran_mod_comp)
c_mod_comp = 'gcc cfunc.c -c -o cfunc.o -O3 -fPIC'
print(c_mod_comp)
system(c_mod_comp)
shared_obj_comp = 'gfortran pygfunc.f90 -c -o pygfunc.o -O3 -fPIC'
print(shared_obj_comp)
system(shared_obj_comp)

ext_modules = [Extension(# module name:
                         'merge',
                         # source file:
                         ['merge.pyx'],
                         # other compile args for gcc
                         extra_compile_args=['-fPIC', '-O3'],
                         # other files to link to
                         extra_link_args=['gfunc.o', 'pygfunc.o','cfunc.o'])]

setup(name = 'merge',
      cmdclass = {'build_ext': build_ext},
      # Needed if building with NumPy.
      # This includes the NumPy headers when compiling.
      include_dirs = [get_include()],
      ext_modules = ext_modules)

table1.py

from finitediff import get_weights
import numpy as np

c = get_weights(np.array([-4.,-3.,-2.,-1.,0,1.,2.,3.,4.]), 0, maxorder=4)
for i in range(len(c[0,:])):
    print('orde (derivative) = ',i)
    print(c[:,i])
    print('')

table2.py

from finitediff import get_weights
import numpy as np

c = get_weights(np.array([(-7./2.),(-5./2.),(-3./2.),(-1./2.),(1./2.),(3./2.),(5./2.),(7./2.)]), 0, maxorder=4)
for i in range(len(c[0,:])):
    print('orde (derivative) = ',i)
    print(c[:,i])
    print('')

table3.py

from finitediff import get_weights
import numpy as np

c = get_weights(np.array([0,1,2,3,4,5,6,7,8]), 0, maxorder=4)
for i in range(len(c[0,:])):
    print('orde (derivative) = ',i)
    print(c[:,i])
    print('')

table4.py

from finitediff import get_weights
import numpy as np

c = get_weights(np.array([(-1./2.),(1./2.),(3./2.),(5./2.),(7./2.),(9./2.),(11./2.),(13./2.),(15./2.)]), 0, maxorder=4)
for i in range(len(c[0,:])):
    print('orde (derivative) = ',i)
    print(c[:,i])
    print('')

test_finitediff.py

from __future__ import print_function, division, absolute_import, unicode_literals

import numpy as np

from finitediff import interpolate_by_finite_diff, derivatives_at_point_by_finite_diff

def test_derivatives_at_point_by_finite_diff():
    out = derivatives_at_point_by_finite_diff(
        np.array([.0, .5, 1.]),
        np.array([.0, .25, 1.]), .5, 2) # y=x**2
    print('')
    print("see derivative approx. in 0.5 for y = x**2 (y' = 2x)")
    print('x:')
    print(np.array([.0, .5, 1.]))
    print('y:')
    print(out)

def test_interpolate_by_finite_diff():
    order = 0
    xarr = np.linspace(-1.5, 1.7, 53)
    yarr = np.exp(xarr)
    xtest = np.linspace(-1.4, 1.6, 57)
    y = interpolate_by_finite_diff(xarr, yarr, xtest)
    yexact = np.exp(xtest)
    for ci in range(y.shape[1]):
        tol = 10**-(13-ci*2)
        print("y = exp(x) so y' = exp(x)")
        print('ytrue:')
        print(yexact)
        print('ycalc:')
        print(y[:,ci])
        print('is approximately similar? ',np.allclose(yexact, y[:,ci]))


if __name__ == '__main__':
    test_interpolate_by_finite_diff()
    test_derivatives_at_point_by_finite_diff()

test_merge.py

from merge import pyadd, pymulti

print(pyadd(2,3))
print(pymulti(2,3))