Numeric Performance in C, C# and Java

By occasion, I find the following draft by Prof. Peter Sestoft:

Numeric Performance in C, C# and Java

He implemented several typical numerical programs in C, C# and Java, and tested their running time.

The conclusion in the paper is that managed code is also performant! Read the paper for the details.

Prof. Sestoft is an expert in compiler design and compiler writing. The draft also contains the generated byte code with his annotations. Enjoy!

Just a note:
Although the C code in the paper is native, it does not take the advantage of the modern CPU design. A fully optimized implementation, which is very hard for normal programmers to write, is faster than the native code. This is why Matlab/Numpy's matrix multiplication is faster than your C code. Matlab uses Intel MKL optimized implementation.

PInvoke native C/C++ libraries in F#

In this blog post, we explore on PInvoke in F#. Before going into PInvoke, an matrix multiplication example is also given to show how to build a DLL for PInvoke, which is important as in most of the cases we don't have pre-built DLLs.

PInvoke stands for Platform Invoke, meaning calling naive libraries. The two major reasons to use PInvoke is 1) native performance, native code is usually faster than .Net managed code, especially for numerical computing routines, and 2) existing stable(well tested) libraries can be used on the fly on .Net.

Of course, it also has some disadvantages, e.g. 1) extra instructions are used for each PInvoke call, this is not important as long as P/Invoke is not used too frequently, 2) calling native code means not safe anymore, e.g. no exception, no memory access checking.

There are several tutorials online for C#: here, here and here. However, these are more focused on existing windows dlls, their motivation for PInvoke is that some functions are not available on .Net, thus calling an existing native DLL, e.g. Win32 API. However, we are talking about performance in this series, PInvoke means that we rewrite the performance critical part in C/C++, or compile an existing library (maybe we also need to make a C wrapper of it in the meantime) and call it using PInvoke. So we should start with the native code part as we don’t have the native dll yet.

Matrix multiplication

Because C and C++ are slightly different, and C++ usually compiles C code well. So we’ll deal with C++ mode in VC++ compiler, i.e. all files have .cpp extensions.

First we write mattrixproduct.cpp:

__declspec(dllexport) int minus(int a, int b)
{
return a - b;
}

extern "C" __declspec(dllexport) int add(int a, int b)
{
return a + b;
}

extern "C" __declspec(dllexport)
void matmul(double *a, int an, int am,
        double *b, int bn, int bm,
        double *c)
{
int i, j, k;

for (i=0; i<an; i++)
    for (j=0; j<bm; j++) {
        double s = 0;
        for (k=0; k<am; k++)
            s += a[i*am + k] * b[k*bm + j];
        c[i*bm+j] = s;
    }
}

__declspec(dllexport) is the header we must add to each function to tell the compiler to generate this function in the DLL. extern "C" is to tell the compiler that this is a C code. We use command line:
cl.exe /LD mattrixproduct.cpp
to get mattrixproduct.dll. We can also create a new Visual C++ project and set the target to .dll.

Here’s the information for this dll: (use command line: link /DUMP /EXPORTS mattrixproduct.dll)

ordinal hint RVA name
1 0 00001000 ?minus@@YAHHH@Z = ?minus@@YAHHH@Z (int __cdecl minus(int,int))
2 1 00001010 add = _add
3 2 00001020 matmul = _matmul

We can see that minus function (without extern “C”) has a non-readable name, while the other two have correct names.

Now, we have the DLL, let’s move on to the PInvoke in F#. Here is the code:

open System
open System.Runtime.InteropServices
open Microsoft.FSharp.NativeInterop
open Microsoft.FSharp.Math

module Native =
[<System.Runtime.InteropServices.DllImport(@"mattrixproduct.dll",
EntryPoint="add")>]
extern int add(int a, int b);

// notice the wired name for minus
[<System.Runtime.InteropServices.DllImport(@"mattrixproduct.dll",
EntryPoint="?minus@@YAHHH@Z")>]
extern int minus(int a, int b);

[<System.Runtime.InteropServices.DllImport(@"mattrixproduct.dll",
EntryPoint="matmul")>]
extern void matmul(double *a, int an, int am, double *b, int bn,
int bm, double *c);


let a = Native.add(10,20)
let b = Native.minus(10,20)
let A = matrix [[1.0;2.0;]; [3.0;4.0;]]
let B = matrix [[3.0;2.0;1.0]; [1.0;1.0;1.0]]
let C = Matrix.zero 2 3
let Ap = PinnedArray2.of_matrix(A)
let Bp = PinnedArray2.of_matrix(B)
let Cp = PinnedArray2.of_matrix(C)
Native.matmul(Ap.Ptr, 2, 2, Bp.Ptr, 2, 3, Cp.Ptr)
Ap.Free()
Bp.Free()
Cp.Free()

printfn "a = %A\nb = %A\nC = %A" a b C

There are several things to mention:

1. For basic type, like int, double, as in the add and minus, they don’t need any special treatment.

2. The hard part is the array. This introduces concept called data marshaling. Because .Net arrays and native arrays are different, so they cannot be directly mapped. This is why we open namespaces for Interop. Another thing is that data marshaling does not support 2 dimensional arrays directly! This is why I used 1-D array in the C implementation. High performance code usually only use 1-D arrays, so this inconvenience does not quite matter to us. To marshal the array, we get a one dimensional pointer for a matrix using PinnedArray2.of_matrix, which is defined in F# PowerPack (in NativeArrayExtensions.fs), and then pass this pointer in the function call.

The performance

Let's now test the performance of the wrapped native matrix multiplication. First write some test on both small matrices and big matrices:

let mm A B =
let Ap = PinnedArray2.of_matrix(A)
let Bp = PinnedArray2.of_matrix(B)
let C = Matrix.zero (A.NumRows) (B.NumCols)
let Cp = PinnedArray2.of_matrix(C)
Native.matmul(Ap.Ptr, A.NumRows, A.NumCols, Bp.Ptr, B.NumRows, B.NumCols, Cp.Ptr)
Ap.Free()
Bp.Free()
Cp.Free()
C

let test() =
// use small matrices
let A = matrix [[1.0;2.0;]; [3.0;4.0;]]
let B = matrix [[3.0;2.0;1.0]; [1.0;1.0;1.0]]
tic()
mm A B |> ignore
toc("my native matrix *")
A * B  |> ignore
toc("F# *")

// use big matrices
let r = new System.Random()
let A1 = Matrix.init (92*112) 280 (fun i j -> r.NextDouble())
let A' = A1.Transpose
tic()
let C1 = mm A' A1
toc("my native matrix *")
let C2 = A' * A1
toc("F# *")

if C1.Equals(C2) then true
else false

The timing functions tic() and toc() are defined in Part IV of the matrix and linear algebra series. Run it in F# interactive:

> test();;
my native matrix * 1.758400 ms
F# * 0.005100 ms
my native matrix * 7207.073800 ms
F# * 16524.670200 ms
val it : bool = true

First the return value is true, indicating that my C implementation is correct. In the small matrix case, the overhead of PInvoke is obvious. In the big matrix case, the performance boost using native code is also obvious.

Note 1: C++ class

So far so good. Because we didn’t touch C++ yet! Classes, templates are more complex.

To work with these. One method is to write C wrappers for C++ classes. This method is commonly used in practice.

We can also directly add __declspec(dllexport) to the class definition. However, creating objects via constructors is not supported anymore. One pattern is to write two static methods: one for creating an object and the other for freeing an object. Member functions are used by providing this pointer explicitly. If you know how to enable object style programming in C, this patter would be familiar to you. Anyway, in this style, the accessing ability is equivalent to C’s. As this pattern requires modifying existing class definitions(although only a little), it is less commonly used, thus I don’t give detailed example here.

Note 2: Pure C

The above example is given in C++ as in most of the cases we are dealing with C++ files, at least C files could be compiled in C++.

For Pure C files (with .c extensions), things are actually easier. Just removing extern "C" would be OK.

Using existing C/C++ code without PInvoke

We can use C++/CLI to compile a C/C++ library without modifying anything into a unsafe managed dll. I will write a separate blog for this later.

Remark

In this tutorial, we know the basics of P/Invoke in F#. We don’t touch two things:

1) Marshaling complex parameters, e.g. different format of strings, non-array pointers, C++ classes, etc.

2) MEMORY! MEMORY! MEMORY! The hard and dangerous part of P/Invoke is memory management. I don’t have enough experience of this part yet. Currently I just follow the style in Math-Provider.

Matrix and linear algebra in F#, Part IV: profile your program, find the bottleneck and speed it up: using matrix multiplication as an example

We will use the face recognition program we developed in Part III to illustrate how to profile programs.

The performance problem of the face recognition program is that it takes about 25 seconds to process 280 faces (each face has 10304=112*92 pixels). I want to find which part takes most of the time. Is it due to the eigen decomposition?

Profile

Profiling means knowing how much time and memory every part of your program take. Good profiling tools could give a detail report on the running behavior of a program. In .Net, there are several good tools for this purpose. Some of them are free.

However, these tools are too big to use in a explorative/interactive programming style. Following Matlab’s tic and toc, I have written them for F#:

let stopWatch = new System.Diagnostics.Stopwatch()

let tic() =
stopWatch.Reset()
stopWatch.Start()

let toc(msg:string) =
stopWatch.Stop()
printfn "%s %f ms" msg stopWatch.Elapsed.TotalMilliseconds
stopWatch.Reset()
stopWatch.Start()

I also use the same style for my C programs:

#include <time.h>
time_t stopWatch;
void tic()
{
stopWatch = clock();
}
void toc(const char *msg)
{
double s = (clock() - stopWatch) * 1000.0 / CLOCKS_PER_SEC;
printf("%s: %.2lf ms\n", msg, s);
stopWatch = clock();
}

This implementation enables us to use a sequence of tocs to profile a code block line by line without the need to write tic every time. However this simple profiler can not be nested as there is only a global timer.

After putting several tocs into my eigCov:

tic()
let colm = colMean faces
toc("mean face")
let A = Matrix.init nrow ncol (fun i j -> faces.[i,j] - colm.[i])
toc("minus mean face")
let L = A.Transpose * A
toc("get L")
let val1, vec1 = eig L 
toc("lapack eig")
let v = val1 * (1.0 / float (ncol - 1))
let u = A * vec1.Transpose
toc("get u")
// normalize eigen vectors
let mutable cnt = 0
for i=0 to ncol-1 do
if abs v.[i] < 1e-8 then
   u.[0..nrow-1,i..i] <- Matrix.create nrow 1 0.0
else
   cnt <- cnt + 1
   let norm = Vector.norm (u.[0..nrow-1,i..i].ToVector())
   u.[0..nrow-1,i..i] <- u.[0..nrow-1,i..i] * ( 1.0 / float norm )
toc("normolize u")

It gets the following output:

mean face 639.891700 ms
minus mean face 195.484600 ms
get L 12885.751500 ms
lapack eig 247.015700 ms
get u 7169.749100 ms
normolize u 623.919700 ms

The bottle neck is in "Get L” and “get u”, both of which are matrix multiplication operations. The main computation part “lapack eig” actually costs very little time (as the 280-by-280 matrix is not big).

The matrix multiplication in “get L” is 280x10304 multiplies 10304x280, which is quite large!

Matrix Multiplication

Matrix multiplication is a great example to show that why we need well tuned libraries to perform numerical computations.

Table: the time cost for self-made(not optimized) implementations

implementations	operator * for matrix type	3 F# loops(use matrix)	3 F# loops(use float Array2D)	3 C loops (native)
time(seconds)	13.07	38.81	13.52	8.10

I used the following 3 loops:

for (i=0; i<N; i++)
for (j=0; j<N; j++) {
   double s = 0;
   for (k=0; k<M; k++)
       s += a[i][k] * b[k][j];
   c[i][j] = s;
}

The two F# programs are similar to the above one in C. We can find that if we use a lot of F# matrix element access operator [,], the performance is very poor. The * operator of matrix type does not have any optimization in it, so it costs the same amount of time as 3-loops (Array2D) does. Native code is faster than the managed ones as it does not do boundary checking for index.

Matlab, R and Numpy/Python all wrap optimized matrix implementations. Let us next see how they perform:

Table: the time cost for optimized implementations

Software	Matlab 2009a	R 2.10.1	Numpy 1.3.0
time	0.25	1.98	0.35

Matlab is the fastest, but it uses 200% CPU. While R and Numpy are single threaded. R takes 2 seconds. Matlab uses Intel MKL, Numpy uses Lapack(maybe optimized version), I don’t know what algebra routines R uses.

We can see the difference between optimized versions and non-optimized ones is very large. A similar report from a .Net math software vender is here. The conclusion is that we should use optimized matrix multiplication procedure when matrices are large.

Add optimized matrix multiplication support to math-provider

The math-provider has made an wrapper for dgemm_, which is a BLAS procedure for computing:

C  :=
alpha*op( A )*op( B ) + beta*C
where op (A) = A or A’. 

This is general case of matrix multiplication. The math-provider also has an wrapper for dgemv_, matrix-vector multiplication.

However, they are not directly callable from Lapack service provider.

add the following code to linear_algebra_service.fs:
let MM a b =

  Service().dgemm_(a,b)

and add the following to linear_algebra.fs:

let MM A B =
if HaveService() then
  LinearAlgebraService.MM A B
else
  A * B

and add the following to linear_algebra.fsi:
/// BLAS matrix multiplication

val MM: A:matrix -> B:matrix -> matrix

Now we define a local name for it:

let mm = LinearAlgebra.MM 

Using BLAS matrix multiplication routine costs about 4 seconds, which is faster than native C’s performance(8 seconds), but worse than Matlab’s or Numpy’s. The reason is that the BLAS(Netlib BLAS 3.1.1) is not optimized for my platform. If I use ATLAS for BLAS and use an Intel Fortran Compiler, the performance will be close to Numpy’s or Matlab’s.

By using BLAS’s matrix multiplication routine for the face recognition, the total running reduces from 25 seconds to 9 seconds.

Attachment

wait.