This report is the continuation of the work done in:

Basic Comparison of Python, Julia, R, Matlab and IDL

Here we:

1. Add new versions of languages
2. Add JAVA
3. Add more test cases.
4. For each language, consistantly use the same method to measure the elapsed time.
5. Provide source codes for all the test cases.
6. Present all the timing results to the fourth digit accuracy (any number less tha 0.0001 is rounded to 0).

While reading this report, be mindful of the following:

• Our intention is not to claim that one language is better than the other.
• In our work, we are often asked to address users’ issues on the computing languages Python, Matlab, IDL, etc. We only have few hours to understand the coding principles of those languages and quickly write codes that resolve users’ issues. We present results in the point of view of a novice programmer.
• If you are an advanced programmer or a language developer and you have results (obtained with optimization techniques) you want to share, feel free to contact us (with a web link) and we wil provide a link to your results here.

All the experiments presented here were done on Intel Xeon Haswell processor node. Each node has 28 cores (2.6 GHz each) and 128 Gb of available memory. We consider the following versions of the languages:

Problem 1: Fibonacci Number

The Fibonacci numbers are the sequence of numbers defined by the linear recurrence equation:

with .

Fibonacci numbers find applications in the fields of economics, computer science, biology, combinatoric, etc.

We implemented both the iterative method and the recursive one, and we record the elepased time for generating the Fibonacci numbers for a given n.

Table 1.1: Elapsed times (in seconds) obtained by computing Fibonacci numbers using then iterative method.

Table 1.2: Elapsed times (in seconds) obtained by computing Fibonacci numbers using the recursive method.

Problem 2: Copy Arrays

This test case is meant to show how fast languages access non-contiguous memory locations.

Consider an arbitrary nxnx3 matrix A. We want to perform the following operations on A:

A(i,j,1) = A(i,j,2)
A(i,j,3) = A(i,j,1)
A(i,j,2) = A(i,j,3)

For instance, in Python the code looks like:

for i in range(n):
     for j in range(n):
          A[i,j,1] = A[i,j,2]
          A[i,j,3] = A[i,j,1]
          A[i,j,2] = A[i,j,3]

The above code segment uses loops. We are also interested on how the same operations are done using vectorization:

A[:,:,1] = A[:,:,2]
A[:,:,3] = A[:,:,1]
A[:,:,2] = A[:,:,3]

The problem allows us to see how each language handles loops and vectorization. We record the elapsed time needed to do the array assignments.

Table 2.1: Elapsed times (in seconds) obtained by copying a matrix using loops.

Table 2.2: Elapsed times (in seconds) obtained by copying a matrix using vectorization.

Problem 3: Matrix Multiplication

We multiply two randomly generated nxn matrices A and B:

C=AxB

This problem shows the importance of taking advantage of built-in libraries available in each language. The elapsed times presented here only measure the times spent on the multiplication (as the size of the matrix varies).

Table 3.1: Elapsed times (in seconds) obtained by multiplying two randomly generated matrices.

Problem 4: Gauss-Legendre Quadrature

The Gauss-Legendre quadrature formulas approximate the integral of a functionby a weighted sum of function-values. When m function-values are used, the formula is exact for polynomials of degree zero through 2m — 1.

Table 4.1: Elapsed times (in seconds) obtained by performing the Gauss-Legendre qudrature.

Problem 5: Numerical Approximation of the 2D Laplace Equation

We find the numerical solution of the 2D Laplace equation:

We use the Jacobi iterative solver. We are interested in fourth-order compact finite difference scheme (Gupta, 1984):

The Jacobi iterative solver stops when the difference of two consecutive approximations falls below 10^{-6}.

Table 5.1: Elapsed times (in seconds) obtained by numerically solving the Posson equation using a Jacobi iterative solver with loops.

Table 5.2: Elapsed times (in seconds) obtained by numerically solving the Posson equation using a Jacobi iterative solver with vectorization.

Problem 6: Belief Propagation

The Belief Propagation can be applied to fields such as speech recognition, computer vision, image processing, medical diagnostics, parity check codes, etc. Its calculations involve a repeated sequence of matrix multiplications, followed by normalization. The Matlab, C and Julia codes are shown in the Justin Domke's weblog (Domke 2012). We report the computing times for various values of the number of iterations (N) when the matrix dimension is 5000x5000.

Table 6.1: Elapsed times (in seconds) obtained by doing the Belief Propagation computations.

Problem 7: Metropolis-Hastings Algorithm

The Metropolis–Hastings (M–H) algorithm is a method for obtaining random samples from a probability distribution. We perform calculations for the implementation of a Metropolis-Hastings algorithm using a two dimensional distribution (Domke 2012). Results are shown when the number of iterations (N) varies.

Table 7.1: Elapsed times (in seconds) obtained by doing the Metropolis algorithm computations.

Problem 8: Manipulation of netCDF Files

We have a set of daily NetCDF files (7305) covering a period of 20 years (1990-2009). The files for a given year are in a sub-directory labeled YYYY (for instance 1990, 1991, 1992, etc.). We want to write a script that  opens each file, reads a three-dimensional variable (longitude/latitude/level), and manipulates it. A pseudo code for the script reads:

Loop over the years
         Obtain the list of NetCDF files
         Loop over the files
                  Read the variable (longitude/latitude/level)
                  Compute the zonal mean average (new array of latitude/level)
                  Extract the column array at latitude 86 degree South
                  Append the column array to a "master" array (or matrix)

The goal here is to be able to do a generate the data to do a contour plot that looks like (obtained with Python):

This is the kind of problems that a typical user we support faces: a collection of thousands of files that needs to be manipulated to extract the desired information. Having tools that allow us to quickly read data from files (in formats such as NetCDF, HDF4, HDF5, grib) is critical for the work we do.

Note that unlike  in Problem 4 of the previous report (where the daily files are in directories associated with months), the daily files to be read in in this case are stored in directories associated with the years. The access to the files is easier in this current problem and we expect the timing numbers to be reduced.

We report in Table 8.1 the elapsed times it took to solve Problem 8 with the various languages.

Table 8.1: Elapsed time (in seconds) obtained by manipulating 7305 NetCDF files on a single processor.

All the above runs were conducted on a node that has 28 cores. Basically, only one core was used. We want to take advantage of all the available cores by spreading the reading of the files and making sure that the data of interest are gathered in the proper order. We use the multi-processing capabilities of the various languages to slightly modify the scripts. For each year, the daily files are read in by different threads (cores).The results are shown in Table 8.2.

Table 8.2: Elapsed time (in seconds) obtained by manipulating 7305 NetCDF files using multiple threading.

Problem 9: Function Evaluations

We  create an array x of length n and loop several times to perform the six operations:

y = sin(x)
x = asin(y)
y = cos(x)
x = acos(y)
y = tan(x)
x = atan(y)

Problem 10: Simple FFT

We create a nxn random complex matrix M and compute the following:

r = fft(M)
r = abs(r)

Problem 11: Square Root of a Matrix

We consider an nxn matrix A with 6s on the diagonal and 1s everywhele else. We are lloking for the matrix B such that BxB = A. We record the time for determining B.

Problem 12: Look and Say Sequence

We write codes to determine the look and say number of order n. Instead of starting with a single digit, we begin with 1223334444.

This test case highlights how languages manipulate strings of arbitrary length.

References

1. Justin Domke, Julia, Matlab and C, September 17, 2012.
2. Michael Hirsch, Speed of Matlab vs. Python Numpy Numba CUDA vs Julia vs IDL, June 2016.
3. Murli M. Gupta, A fourth Order poisson solver, Journal of Computational Physics, 55(1):166-172, 1984.
4. Jean Francois Puget, A Speed Comparison Of C, Julia, Python, Numba, and Cython on LU Factorization, January 2016.
5. Alex Rogozhnikov, Log-likelihood benchmark, September 2015.
6. Sebastian Raschka, Numeric matrix manipulation - The cheat sheet for MATLAB, Python Nympy, R and Julia, June 2014.
7. Yousef Saad, Iterative Methods for Sparse Linear Systems (2 ed.), SIAM, ISBN 0898715342, 200366

Source Files

All the source files for the problems presented here are in the attached file: sourceFiles2018.tar.gz

If you have a comment/suggestion/question, contact Jules Kouatchou (Jules.Kouatchou@nasa.gov)gcc