Skip to content

Instantly share code, notes, and snippets.

@yudhastyawan

yudhastyawan/cfunc.c

Last active April 25, 2021 04:55
Show Gist options
  • Save yudhastyawan/fcc18cde076869905d54ba9b8b9aed71 to your computer and use it in GitHub Desktop.
Save yudhastyawan/fcc18cde076869905d54ba9b8b9aed71 to your computer and use it in GitHub Desktop.
Download ZIP
homework 3: embedded fortran and C in Python
Raw
cfunc.c
#include <stdio.h>
#include "cfunc.h"
void c_add(double* a, double* b, double* c)
{
*c = *a + *b;
}
Raw
cfunc.h
void c_add(double* a, double* b, double* c);
Raw
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
Raw
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
Raw
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
Raw
pygfunc.h
extern void c_multi(double* a, double* b, double* c);
Raw
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)
Raw
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('')
Raw
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('')
Raw
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('')
Raw
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('')
Raw
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()
Raw
test_merge.py
from merge import pyadd, pymulti
print(pyadd(2,3))
print(pymulti(2,3))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment