Writing kernel modules in Haskell [LWN.net]

Writing kernel modules in Haskell

Posted Sep 15, 2009 10:35 UTC (Tue) by sagrailo (guest, #48054) [Link]

Yeah, right... I just quickly wrote this code (to multiply two matrices):

#include <algorithm>
#include <cassert>
#include <cstdlib>
#include <ctime>
#include <iostream>

extern "C" {
#include <cblas.h>
}

#include <tnt.h>

#define N 1024

int
main(void)
{
    float* A = new float[N * N];
    float* B = new float[N * N];
    float* C = new float[N * N];
    std::fill(C, C + N * N, 0);

    TNT::Array2D<float< ATNT(N, N);
    TNT::Array2D<float< BTNT(N, N);
    TNT::Array2D<float< CTNT(N, N);

    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j) {
            A[i * N + j] = ATNT[i][j] = 1;
            B[i * N + j] = BTNT[i][j] = 1;
        }

    clock_t start, end;

    start = clock();
    cblas_sgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, N, N, N, 1, A, N, B,
                N, 0, C, N);
    end = clock();
    std::cerr
        << "ATLAS: "
        << (float) (end - start) / CLOCKS_PER_SEC
        << "s"
        << std::endl;

    start = clock();
    CTNT = TNT::matmult(ATNT, BTNT);
    end = clock();
    std::cerr
        << "TNT:   "
        << (float) (end - start) / CLOCKS_PER_SEC
        << "s"
        << std::endl;

    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j) {
            assert(C[i * N + j] == N);
	    assert(CTNT[i][j] == N);
        }

    delete [] A;
    delete [] B;
    delete [] C;
}

Complied with (I've put TNT into /tmp, and I have ATLAS installed in /usr/local):

g++ -O3 -o foo -I/tmp/tnt foo.cpp -lcblas -latlas

The output, on my old ThinkPad T61p:

ATLAS: 0.22s
TNT:   17.64s

So - C++ almost two orders of magnitude slower (and that's typical, in my experience).

Writing kernel modules in Haskell

Posted Sep 15, 2009 12:01 UTC (Tue) by jbh (subscriber, #494) [Link]

Interesting! I've wanted to look into Eigen as well, which is another
template linear algebra library (http://eigen.tuxfamily.org). It seems to
hold up pretty well against Atlas (on Thinkpad T60, compiled with -msse2
flag):

ATLAS: 0.43s
TNT: 32.25s
Eigen: 0.76s

The major additions to your test program are (in different places):

#include <Eigen/Core>
--
MatrixXf Aeig(N,N), Beig(N,N), Ceig(N,N);
Aeig.setOnes(); // or Aeig = MatrixXf::Ones(N,N);
Beig.setOnes();
Ceig.setZero();
--
Ceig = Aeig * Beig;
--
assert(Ceig(i,j) == N);

Writing kernel modules in Haskell

Posted Sep 15, 2009 12:50 UTC (Tue) by bjacob (subscriber, #58566) [Link]

Hm, small world, I'm actually one of the main developers of Eigen. I didn't want to make a plug, but now that you've mentioned it...

I'm surprised by your result, because in our experience, we are consistently faster than ATLAS for matrix products. Already with Eigen 2.0, and even more so in the development branch. In fact, on single-CPU systems, we have almost the performance of vendor-tuned BLAS.

development branch:
http://eigen.tuxfamily.org/index.php?title=Benchmark

2.0:
http://eigen.tuxfamily.org/index.php?title=Benchmark-Augu...

It's good that you passed -msse2, indeed.

Some things to consider:

* did you allow ATLAS to use more than one thread? That would probably not be applicable in the setting of a kernel module! (As far as I understand) so perhaps it's more fair to require ATLAS to use only one thread.

* did you tune ATLAS for your hardware, especially CPU cache size? You can do the same for Eigen, then. For example with the development version it's
EIGEN_TUNE_FOR_CPU_CACHE_SIZE, see in the file Macros.h.

Writing kernel modules in Haskell

Posted Sep 15, 2009 13:09 UTC (Tue) by jbh (subscriber, #494) [Link]

In short, no tuning at all, neither ATLAS or Eigen. Standard ubuntu (jaunty)
packages for everything. This was on my laptop, and I only bother to tune
for production runs; those are not on my laptop :)

Eigen already looks very impressive, for a relatively young project in a
problem space that's dominated by dinosaurs --- I've meant to look into it
for a while, but not enough time and all that, and last time I checked it
didn't have sparse matrix support (I see it does now, although
experimental).

Writing kernel modules in Haskell

Posted Sep 15, 2009 13:29 UTC (Tue) by bjacob (subscriber, #58566) [Link]

OK, my best guess then is that ATLAS by default uses more than one thread...
anyway thanks for the kind words. Yes, there are a lot of dinosaurs in this area and we're sort of the small mammals in comparison.

(If it wasn't obvious, the largest advantage of a c++ template library like Eigen over a large C or Fortran binary library like a BLAS, is that it spontaneously generates only the code that's needed at build-time, and after that you can forget about it, so your kernel module is only a few dozen kilobytes self-contained; by comparison vendor-tuned BLAS take dozens of megabytes).

Writing kernel modules in Haskell

Posted Sep 29, 2009 23:14 UTC (Tue) by Auders (subscriber, #53318) [Link]

I did the same test with standard jaunty packages and the same compilation parameters and for me Eigen outperformed ATLAS, as you expected:
ATLAS: 0.69s
TNT: 21.93s
Eigen: 0.56s

Matrix Libraries

Posted Sep 16, 2009 20:50 UTC (Wed) by cry_regarder (subscriber, #50545) [Link]

Of course ublas has bindings for atlas so you can use this approach:

...
#include <boost/numeric/ublas/matrix.hpp>
#include <boost/numeric/bindings/atlas/cblas.hpp>
#include <boost/numeric/bindings/traits/ublas_matrix.hpp>
...
namespace ublas = boost::numeric::ublas;
namespace atlas = boost::numeric::bindings::atlas;
...
atlas::gemm(CblasNoTrans, CblasNoTrans, 1.0,
Aublas, Bublas, 0.0, Cublas);
...

and this compile:
g++ -O3 -o mat mat.c++ -I/tmp/boost-numeric-bindings -L /usr/lib64/atlas -latlas -lcblas -DNDEBUG

to get run times like:

ATLAS: 0.16s
uBLAS: 0.16s

...

Cry