PROGRAMMING LANGUAGE INTERSTICES ROHIT GOSWAMI Created: 2021-09-09 Thu 12:52

BRIEF INTRODUCTION 3

HELLO! Find me here: Who? Rohit Goswami MInstP Doctoral Researcher, University of Iceland, Faculty of Physical Sciences https://rgoswami.me 4

LOGISTICS All contents are Slides are in presentations/lesHouches21 hosted on GitHub Questions are welcome anytime 5

PROGRAMMING LANGUAGES 7

MOTIVATION –> Through the magic of automated coding and standards “If a program or package (the words are used interchangeably) is to have a long life and to be of wide application in its field, it is essential for it to be easily moved from one machine to another. It used to be common to dismiss such movement with the statement, ‘There is no such thing as a machine-independent program.’ Nonetheless, a great many packages do now move from one machine to another”[lyonUsingAnsFortran1980] 8

LANGUAGE STANDARDS “The standard is the contract between the compiler writer and the application developer.”[clermanModernFortranStyle2012] Sadly interpreters are not really Code -> Memory Python itself is an interpreter which is commonly implemented in C 9

CHANGING STANDARDS character(10) BLAH*8 character*8 :: BLAH_ONE(10) character(8) :: BLAH_ONE(10) #!/usr/bin/env python print("Hello World") print "Hello World" 10

F77 ∉ F90 ALWAYS

LFORTRAN 13

INTRODUCTION Language Files Lines Code Comments Blanks C 3 1023 694 191 138 C Header 57 14237 11416 1041 1780 CMake 11 430 361 16 53 C++ 54 30745 26911 1560 2274 C++ Header 1 8938 8297 348 293 FORTRAN 20 1738 1344 174 220 Python 2 224 191 4 29 Total 148 57335 49214 3334 4787 14

STRUCTURE

FEATURES Runtime Libraries Pure Fortran | Impure ASR Guarantees compilation –> Wrappers Parser Currently almost all of F2018; including F77 LLVM Canonical backend, includes compile time evaluated expressions Jupyter Native execution as a kernel 16

PRACTICAL MATTERS 18

NUMERICAL METHODS 101 Numerical stability issues Floating Point Numbers IEEE 754 Convergence can be tied to machine epsilon…. Or a magic number Or Chemical/Spectroscopic accuracy 𝐴𝑥 = 𝑏 𝑥 = 𝑏 𝐴 −1 [goldbergWhatEveryComputer1991] (1𝑒 − 5, 1𝑒 − 16) program main print*, tiny(0.d0) end program 2.2250738585072014e-308 19

YAWN “I don’t care because I use libraries. ” “I don’t care because I write my own algorithms.” 20

IMPLEMENTING SINE real(dp) function dsin(x) result(r) real(dp), parameter :: pi = 3.1415926535897932384626433832795_dp real(dp), intent(in) :: x ! Assume folded to [-2𝜋 ,2𝜋 ] real(dp) :: y, xnew, twoPi, invTwoPi if (abs(x) < pi/2) then r = kernel_dsin(x) else ! fold to pi/2 from https://realtimecollisiondetection.net/blog/?p= y = modulo(xnew, 2*pi) y = min(y, pi - y) y = max(y, -pi - y) y = min(y, pi - y) r = kernel_dsin(y) end if end function

FOLDING I Symmetric about –> 𝜋/2 𝑥 − 𝜋/2 𝑥 = −𝑎𝑏𝑠(𝑥)

FOLDING II -> 𝑥 = −𝑎𝑏𝑠(𝑥– 𝜋/2) + 𝜋/2 𝑥 = 𝑚𝑖𝑛(𝑥, 𝜋 − 𝑥) 23

FOLDING III Repeat rightward to get Complete with 𝑥 = 𝑚𝑎𝑥(𝑥, −𝜋 − 𝑥) 𝑥 = 𝑚𝑖𝑛(𝑥, 𝜋 − 𝑥)

KERNEL SINE [ChevillardJoldesLauter2010] nix-shell -p sollya rlwrap rlwrap -A sollya prec=120; f=sin(x); d=[-pi/2;pi/2]; # Use min/max poly p=remez(f,16,d, default,1e-16,1e-30); p; # Returns coefficients # Zero out even terms ! Accurate on [-pi/2,pi/2] to about 1e elemental real(dp) function kernel_dsi real(dp), intent(in) :: x real(dp), parameter :: S1 = 0.99999999 real(dp), parameter :: S2 = -0.1666666 real(dp), parameter :: S3 = 8.33333333 real(dp), parameter :: S4 = -1.9841269 real(dp), parameter :: S5 = 2.75573155 real(dp), parameter :: S6 = -2.5051823 real(dp), parameter :: S7 = 1.60465859 real(dp), parameter :: S8 = -7.3572396 real(dp) :: z z = x*x res = x * (S1+z*(S2+z*(S3+z*(S4+z*(S5+ (S6+z*(S7+z*S8))))))) end function 25

AUXILIARY FUNCTIONS elemental real(dp) function dmodulo(x, y) result(r) real(dp), intent(in) :: x, y r = x-floor(x/y)*y end function elemental integer function dfloor(x) result(r) real(dp), intent(in) :: x if (x >= 0) then r = x else r = x-1 end if end function elemental real(dp) function dabs(x) result(r) real(dp), intent(in) :: x if (x >= 0) then r = x else r = -x end if end function 26

ACCURACY With an evenly spaced grid of 0.01: over kernel range 1𝐸 − 16 −𝜋/2, 𝜋/2 over 1𝐸 − 16 [−20, 20] from 1𝐸 − 12 1𝐸5 from 1𝐸 − 7 1𝐸10 27

POST MORTEM The problem turned out to be modulo IEEE defines a remainder for trignometric functions ieee754_rem_pio2 Do not naively assume “handcrafted” algorithms are IEEE 754 accurate 28

FORTRAN, C, PYTHON F2003 Introduced the ISO_C_BINDING (no term) Described in great detail on : fortran90.org 30

F2PY 32

HISTORY Developed by Pearu Peterson July 9, 1999 f2py.py –> Fortran to Python Interface Generator (FPIG) January 22, 2000 f2py2e –> Fortran to Python Interface Generator, 2nd edition. July 19, 2007 numpy.f2py –> f2py2e moved to NumPy project. This is current stable code of f2py. Used extensively for F77 , :) [petersonF2PYToolConnecting2009] NumPy [waltNumPyArrayStructure2011] SciPy [virtanenSciPyFundamentalAlgorithms2020] MsSpec [sebilleauMsSpec1MultipleScattering2011] 33

DESIGN A best effort wrapper Specifications via .pyf or inline comments Not a compiler Can rewrite code :) 34

ROADMAP Moves about as fast as spectroscopy codes :) Also picking up speed in 2021 Implementing newer standards (90, 95, 2003, 2008, 2018, 2020Y) Automating guarantees Documentation 35

RELEVANCE Writing efficient wrappers without being a language lawyer 36

BONUS: READING CODE 38

READING CODE I main: push rbp mov rbp, rsp mov DWORD PTR [rbp-4], 3 mov eax, 0 pop rbp ret __static_initialization_ and_destruction_0(int, int): push rbp mov rbp, rsp sub rsp, 16 mov DWORD PTR [rbp-4], edi mov DWORD PTR [rbp-8], esi cmp DWORD PTR [rbp-4], 1 jne .L5 cmp DWORD PTR [rbp-8], 65535 jne .L5 mov edi, OFFSET FLAT:_ZStL8 __ioinit call std::ios_base::Init::Init() [complete object constructor] mov edx, OFFSET FLAT:__dso_handle mov esi, OFFSET FLAT:_ZStL8__ioini mov edi, OFFSET FLAT:_ZNSt8ios_bas call __cxa_atexit .L5: nop leave ret _GLOBAL__sub_I_main: push rbp mov rbp, rsp mov esi, 65535 mov edi, 1 call __static_initialization_ and_destruction_0(int, int) pop rbp ret But who writes assembly anyway? 39

READING CODE II int main () { int D.48918; { int a; a = 3; D.48918 = 0; return D.48918; } D.48918 = 0; return D.48918; } void _GLOBAL__sub_I_main.cpp () { __static_initialization_ and_destruction_0 (1, 65535); } void __static_initialization_ and_destruction_0 (int __initialize_p, int __priority) { if (__initialize_p == 1) goto ; : if (__priority == 65535) goto ; : std::ios_base::Init::Init (&__ioinit __cxxabiv1::__cxa_atexit (__dt_comp &__ioinit, &__dso_han goto ; : : goto ; : : } GIMPLE is an internal gcc representation… 40

READING CODE III #include int main() { int a=3; return 0; } Better for most people, still a bit lacking for novices Assigning an integer Produces a file binary which can be run as: Output There is no output, but an assignment of an integer with value 3 takes place g++ main.cpp -o file ./file What about different languages? 41

READING CODE IV Maybe gcc is just an ugly compiler… program main integer :: x = 3 + 6 print *, x end program lfortran has a nicer intermediate structure conda create -n lf conda activate lf conda install lfortran \ -c conda-forge lfortran --show-asr consint.f90 42

CONCLUSIONS 44

OMITTED TOPICS Web development and design Including frameworks and UX Continuous integration How to ensure documentation is coupled to working code Benchmarking Demonstrating code superiority Code Review Practices Scrum and teamwork HPC and Parallelism Efficient data usage and algorithms Build systems Portability Science Also algorithms Usage examples For f2py and lfortran 45

KEY TAKEAWAYS Document at every level Write modern Fortran Use ISO_C_BINDINGS Use IEEE standards When using F77 Use f2py Don’t be too clever Don’t be too accurate Ensure style guides exist Benchmark Lint automatically Use modern tools , fpm stdlib LFortran Ask for help -> Fortran Lang 46

THE END 48

ACKNOWLEDGEMENTS as my supervisor, as my co-supervisor, and my committee member at Los Alamos National Laboratory ( , and ) Family, pets, Groupmembers, audience Prof. Hannes Jónsson Prof. Birgir Hrafnkelsson Dr. Elvar Jonsson Dr. Ondřej Čertík Quansight Labs Dr. Ralf Gommers Dr. Melissa Weber Mendonça Dr. Pearu Peterson 49

BIBLIOGRAPHY Chevillard, Joldeş & Lauter, Sollya: An Environment for the Development of Numerical Codes, 28-31, in in: Mathematical Software - ICMS 2010, edited by Fukuda, van der Hoeven, Joswig & Takayama, Springer Clerman & Spector, Modern Fortran: Style and Usage, Cambridge University Press . Goldberg, What Every Computer Scientist Should Know about Floating-Point Arithmetic, ACM Computing Surveys, 23(1), 5-48 . . . Lyon, Using Ans Fortran, National Bureau of Standards . Peterson, F2PY: A Tool for Connecting Fortran and Python Programs, International Journal of Computational Science and Engineering, 4(4), 296 . . . Sébilleau, Natoli, Gavaza, Zhao, Da Pieve & Hatada, MsSpec-1.0: A Multiple Scattering Package for Electron Spectroscopies in Material Science, Computer Physics Communications, 182(12), 2567-2579 . . . Virtanen, Gommers, Oliphant, Haberland, Reddy, Cournapeau, Burovski, Peterson, Weckesser, Bright, van der Walt, Brett, Wilson, Millman, Mayorov, Nelson, Jones, Kern, Larson, Carey, Polat, Feng, Moore, VanderPlas, Laxalde, Perktold, Cimrman, Henriksen, Quintero, Harris, Archibald, Ribeiro, Pedregosa & van Mulbregt, SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python, Nature Methods, 17(3), 261-272 . . . van der Walt, Colbert & Varoquaux, The NumPy Array: A Structure for Efficient Numerical Computation, Computing in Science Engineering, 13(2), 22-30 . . [ChevillardJoldesLauter2010] [clermanModernFortranStyle2012] [goldbergWhatEveryComputer1991] link doi [lyonUsingAnsFortran1980] [petersonF2PYToolConnecting2009] link doi [sebilleauMsSpec1MultipleScattering2011] link doi [virtanenSciPyFundamentalAlgorithms2020] link doi [waltNumPyArrayStructure2011] doi 50

THANKS! 51