I'm trying to add Eigen into kind of quick-and-dirty SGEMM profiler of mine. Basically, I'm allocating three n*n arrays, am then initializing first two of them, and finally am calling corresponding matrix mutliplication routines. Here is one I wrote for using Eigen:

Code: Select all

float
sgemmEigen(const int n, const float *a, const float *b, float *c)
{
    std::clock_t start = std::clock();

    Eigen::Map < Eigen::MatrixXf > am(a, n, n);
    Eigen::Map < Eigen::MatrixXf > bm(b, n, n);
    Eigen::Map < Eigen::MatrixXf > cm(c, n, n);
    cm = am * bm;

    std::clock_t stop = std::clock();
    return (float) (stop - start) / CLOCKS_PER_SEC;
}


However, results returned by this routine seem to be wrong; so I guess I'm not using Map<> properly - any hint here?

Just in case, here is this code put together with the test code I'm using, so that this one is complete program to try:

Code: Select all

#include <cassert>
#include <ctime>
#include <iostream>

#include <Eigen/Core>

#define N 768
#define EPSILON 1e-4

int
main(void)
{
        float          *a = new float[N * N];

        for (int k = 0, i = 0; i < N; ++i)
                for (int j = 0; j < N; ++j, ++k)
                        a[k] = (float) (i + 1) / N;

        float          *b = new float[N * N];
        assert(b != NULL);
        for (int k = 0, i = 0; i < N; ++i)
                for (int j = 0; j < N; ++j, ++k)
                        b[k] = (float) (j + 1) / N;

        float          *c = new float[N * N];
        assert(c != NULL);

        std::clock_t start = std::clock();

        Eigen::Map < Eigen::MatrixXf > am(a, N, N);
        Eigen::Map < Eigen::MatrixXf > bm(b, N, N);
        Eigen::Map < Eigen::MatrixXf > cm(c, N, N);
        cm = am * bm;

        std::clock_t stop = std::clock();
        float elapsed = (float) (stop - start) / CLOCKS_PER_SEC;
        std::cout << N * N * 2 * N / elapsed / 1e9 << " GFLOPS/s\n";

        for (int k = 0, i = 0; i < N; ++i)
                for (int j = 0; j < N; ++j, ++k)
                        assert(fabs(c[k] - (float) (i + 1) * (j + 1) / N) /
                               ((float) (i + 1) * (j + 1) / N) < EPSILON);

        delete [] a;
        delete [] b;
        delete [] c;

        return 0;
}


I'm on 64-bit Slackware 13.0, Eigen is last update available for my distro (and that would be 2.0.10), gcc is 4.4.3.

I confirm that you're using Map correctly. Can't say much more, sorry

bjacob wrote:I confirm that you're using Map correctly. Can't say much more, sorry

I tried with simpler code example:

Code: Select all

#include <iostream>
#include <Eigen/Core>

int
main(void)
{
        float          *a = new float[4];
        a[0] = 1, a[1] = 1, a[2] = 2, a[3] = 2;

        float          *b = new float[4];
        b[0] = 1, b[1] = 2, b[2] = 1, b[3] = 2;

        float          *c = new float[4];

        Eigen::Map < Eigen::MatrixXf > am(a, 2, 2);
        Eigen::Map < Eigen::MatrixXf > bm(b, 2, 2);
        Eigen::Map < Eigen::MatrixXf > cm(c, 2, 2);
        cm = am * bm;

        std::cout << c[0] << " " << c[1] << std::endl;
        std::cout << c[2] << " " << c[3] << std::endl;

        delete [] a;
        delete [] b;
        delete [] c;

        return 0;
}


So I'm multiplying matrices (in Matlab notation) [1 1; 2 2] and [1 2; 1 2], and the result should be matrix [2 4; 4 8], but the code above is producing result [5 5; 5 5]. So I'd suspect Map<> is supposing column-major order of underlying data - is that correct?

Map<MatrixXf> interpretes the array as the storage of a MatrixXf, and that is column-major.

If you want to map an array as a row-major matrix, you can do:

Code: Select all

typedef Matrix<float,Dynamic,Dynamic,RowMajor> RMatrixXf;

and then use

Code: Select all

Map<RMatrixXf>

If you use Map on row-major matrix types, you are strongly encouraged to upgrade to the latest revision of the 2.0 branch (to become 2.0.12) as grave bugs exactly in this area were just fixed:
http://eigen.tuxfamily.org/index.php?ti ... e#Download
http://bitbucket.org/eigen/eigen/get/2.0.tar.bz2

You can even tell Matrix to default to row-major order, hence MatrixXf be row-major, by defining the EIGEN_DEFAULT_TO_ROW_MAJOR token.

Finally, since you're measuring performance, note that the development branch would be much faster there, for 2 reasons,
1) matrix-matrix product is faster there,
2) it allows you to map a 16-byte aligned array and tell Eigen that it is aligned, which further helps with SSE vectorization.

Thanks bjacob, that worked for me.

Performance related question: for example on Core 2 Duo machine - is Eigen able to utilize both SSE units? With my quick measurings, as above, Eigen performance is approximately half of GFLOPS available for given machine...

Eigen doesn't do any parallelization at the moment. So each thread that uses Eigen, uses only one of your 2 cores. Eigen doesn't spawn threads by itself.

On the flip side, Eigen is completely stateless so in particular it is completely thread-safe, so if you have 2 independent computations to run, you can run them simultaneously in 2 separate threads.

Hi,

Using the same idea I tried this:

Code: Select all

typedef Array<float, Dynamic, Dynamic, RowMajor> RArrayXXf;
typedef Map<RArrayXXf> RExtMat;

Now, my function goes like that:

Code: Select all

void InitCentroidsPlusPlus(float *mClusterCentroids, float *mDataSamples, int numCentroids, int numSamples, int seedNumber) {
   RExtMat X(mDataSamples, numSamples, NUM_FEATURES);
   RExtMat Mu(mDataSamples, numSamples, NUM_FEATURES);

   SomeCode...
}

Yet using Eiegn 3.3.2 and Visual Studio it creates error.

Any idea?

What are the errors? This should work. (unless mClusterCentroids turns out to be const)

It just crashes.

Nothing is const (I write my Code in C Style, so no const).

ah, then check the backtrace, and also make sure to try it in debug mode to get Eigen's assertions. There is nothing wrong in your code snippet.