Eigen min/max/cast

#include <iostream>
#include <Eigen/Core>
#include <boost/timer.hpp>

using namespace std;
using namespace Eigen;

typedef unsigned short T;
typedef Matrix<T, 3, 1> Vector3T;

const float offset = 3.0f;
const float clamp = 16.0f;

// 297 ms
void eigen1(const Matrix3f & mat, T* in, size_t size, T* out) {
    for (size_t i = 0; i < size*4; i += 4) {
        Map<Vector3T> src(in + i);
        Map<Vector3T> dest(out + i);
        dest = ((mat * src.cast<float>()).array() + offset).max(0.0).min(clamp).cast<T>();
    }
}

// 255 ms
void eigen2(const Matrix3f & mat, T* in, size_t size, T* out) {
    for (size_t i = 0; i < size*4; i += 4) {
        Vector3f src(in[i + 0], in[i + 1], in[i + 2]);
        Map<Vector3T> dest(out + i);
        dest = ((mat * src).array() + offset).max(0.0).min(clamp).cast<T>();
    }
}

// 249 ms
void eigen3(const Matrix3f & mat, T* in, size_t size, T* out) {
    for (size_t i = 0; i < size*4; i += 4) {
        Vector3f src(in[i + 0], in[i + 1], in[i + 2]);
        Vector3f result = ((mat * src).array() + offset).max(0.0).min(clamp);
        out[i + 0] = (T) result(0);
        out[i + 1] = (T) result(1);
        out[i + 2] = (T) result(2);
    }
}

// 152 ms
void eigen4(const Matrix3f & mat, T* in, size_t size, T* out) {
    for (size_t i = 0; i < size*4; i += 4) {
        Vector3f src(in[i + 0], in[i + 1], in[i + 2]);
        Vector3f result = mat * src;
        out[i + 0] = (T) max(0.0f, min(clamp, result(0) + offset));
        out[i + 1] = (T) max(0.0f, min(clamp, result(1) + offset));
        out[i + 2] = (T) max(0.0f, min(clamp, result(2) + offset));
    }
}

// 163 ms
void manual(const Matrix3f & mat, T* in, size_t size, T* out) {
    for (size_t i = 0; i < size*4; i += 4) {
        float x = in[i + 0];
        float y = in[i + 0];
        float z = in[i + 0];
        out[i + 0] = (T) max(0.0f, min(clamp, x*mat(0,0) + y*mat(0,1) + z*mat(0,2) + offset));
        out[i + 0] = (T) max(0.0f, min(clamp, x*mat(0,0) + y*mat(0,1) + z*mat(0,2) + offset));
        out[i + 0] = (T) max(0.0f, min(clamp, x*mat(0,0) + y*mat(0,1) + z*mat(0,2) + offset));
    }
}

int main(int argc, char** argv) {
    Matrix3f mat;
    mat << 10, -1, -2,
           -1, 10, -2,
           -1, -2, 10;

    const size_t size = 10000000;
    T* in = new T[size * 4];
    T* out = new T[size * 4];
    for (int i = 0; i < size * 4; ++i) {
        in[i] = i % std::numeric_limits<T>::max();
        out[i] = 0;
    }

    boost::timer timer;

    if (argc == 2 && strcmp(argv[1], "eigen1") == 0)
    {
        eigen1(mat, in, size, out);
    }
    else if (argc == 2 && strcmp(argv[1], "eigen2") == 0)
    {
        eigen2(mat, in, size, out);
    }
    else if (argc == 2 && strcmp(argv[1], "eigen3") == 0)
    {
        eigen3(mat, in, size, out);
    }
    else if (argc == 2 && strcmp(argv[1], "eigen4") == 0)
    {
        eigen4(mat, in, size, out);
    }
    else if (argc == 2 && strcmp(argv[1], "manual") == 0)
    {
        manual(mat, in, size, out);
    }
    else
    {
        cout << "Usage: " << argv[0] << " eigen1|eigen2|eigen3|eigen4|manual" << endl;
        exit(1);
    }

    double elapsed = timer.elapsed();
    cout << elapsed << endl;
}